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Il 25 giugno 1998, a circa due anni dalla scomparsa, si è tenuta la 
Giornata in ricordo di Placido Cicala organizzata dai Dipartimenti di 
Ingegneria Strutturale e di Ingegneria Aeronautica e Spaziale del Politecnico 
di Torino, col patrocinio dell’Accademia Nazionale dei Lincei, 
dell’Accademia delle Scienze di Torino e della Città di Torino. Hanno 
inoltre concesso un contributo per la realizzazione dell’iniziativa il 
Politecnico di Torino ed il Consiglio Nazionale delle Ricerche. 

La manifestazione si è svolta nella sala del Consiglio di Facoltà del 
Politecnico. 

La Giornata è stata aperta dal Prof. Ugo ROSSETTI, decano dei docenti 
del Dipartimento di Ingegneria Strutturale, che ha successivamente svolto la 
funzione di moderatore. 

È stata anzitutto data lettura del telegramma di adesione dal Prof. 
Edoardo VESENTINI, Presidente dell’Accademia Nazionale dei Lincei, cui 
ha fatto seguito il Prof. Elio CASETTA, Presidente dell’Accademia delle 
Scienze di Torino che ha ricordato Placido CICALA, Socio corrispondente a 
partire dal 1940 e Socio nazionale dal 1953. 

Per il Politecnico, in assenza del Rettore ZICH, ha portato l’adesione il 
Prof. Pietro APPENDINO che ha poi sviluppato il suo intervento in qualità di 
Preside della 18 Facoltà di Ingegneria. 

Hanno successivamente preso la parola il Prof. Piero MARRO, Direttore 
del Dipartimento di Ingegneria Strutturale che ha illustrato l’attività del 
Prof. Cicala nel Dipartimento a partire dal 1956 ed il Prof. Sergio CHIESA, 
Direttore del Dipartimento di Ingegneria Aerospaziale, che ha ricordato 
l’attività del prof. Cicala dal 1934 al 1956 presso l’allora Istituto di 
Costruzioni Aeronautiche. 

Sono successivamente seguiti gli interventi del Prof. Silvio NOCILLA del 
Politecnico di Torino e del Prof. Edoardo BENVENUTO, Preside della Facoltà 
di Architettura dell’Università di Genova, che ha portato il saluto dei 
colleghi docenti di Scienza delle Costruzioni delle altre sedi italiane. 

Il Prof. Ernesto VALLERANI ha poi recato l’adesione dell’ Associazione 
Italiana di Aeronautica ed Astronautica. 


* . . 
Questo volume è stato pubblicato con fondi provenienti da un contributo erogato dal Consiglio 
Nazionale delle Ricerche per l’organizzazione della Giornata. 
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Si sono poi susseguiti i ricordi di tre allievi del Prof. Cicala: la Prof. 
Angia SASSI e il Prof. Franco ALGOSTINO che hanno ricordato il Maestro e 
l’Amico e l’atmosfera serena di collaborazione che Egli aveva saputo creare 
sia nella didattica che nella ricerca, e l’Ing. Dante BIGI che fu il primo 
dottore di ricerca che conseguì il titolo sotto la guida del Prof. Cicala. 

Hanno inoltre inviato la loro adesione il Sindaco di Torino Prof. 
Valentino CASTELLANI ed il Prof. Dionigi GALLETTO dell’Università di 
Torino. 

Al termine della seduta iniziale consacrata ai ricordi di colleghi e di 
allievi, ha avuto inizio la seduta dedicata alle memorie scientifiche in 
programma, il cui testo è riportato nelle pagine seguenti. 





Introduction 


On 25 June, almost two years after his death, the P/acido Cicala 
Memorial Day was held. This meeting was organised by the Department of 
Structural Engineering and the Department of Aeronauties and Space 
Engineering of the Politecnico di Torino, under the Patronage of the 
Accademia Nazionale dei Lincei, the Accademia delle Scienze di Torino 
and the City of Turin. The Politecnico di Torino and the CNR (National 
Research Council) also kindly contributed to this initiative. 

The celebration was held in the Faculty Council Hall of the Politecnico 
di Torino. 

The Memorial Day was opened by Prof. Ugo ROSSETTI, the senior 
member of the academic staff of the Department of Structural Engineering, 
who then acted as chairman. 

First of all a telegram from Prof. Edoardo VESENTINI, President of the 
Accademia Nazionale dei Lincei, was read. There followed an address by 
Prof. Elio CASETTA, President of the Accademia delle Scienze di Torino, 
who recalled Placido Cicala as correspondent Member from 1940 and 
National Member from 1953. 

Prof. Pietro APPENDINO, as Dean of the 1% Faculty of Engineering, 
represented the Politecnico, in the absence of the Rector, Prof. ZICH. 

An address by Prof. Piero MARRO, Director of the Department of 
Structural Engineering, then followed, in which he explained Prof. Cicala’s 
activities in the Department starting from 1956. Prof. Sergio CHIESA, 
Director of the Department of Aerospace Engineering then took the floor. 
He recalled the activity of Prof. Cicala in the Aeronautical Construction 
Institute, from 1934 to 1956. 

Prof. Silvio NOCILLA, from the Politecnico di Torino then spoke and he 
was followed by Prof. Edoardo BENVENUTO, Dean of the Faculty of 
Architecture in Genoa, who brought greetings from the academic colleagues 
of Scienza delle Costruzioni of the other Italian Universities. 

Prof. Ernesto VALLERANI then spoke on behalf of the Associazione 
Italiana di Aeronautica ed Astronautica. 

There then followed the recollections of three of Prof. Cicala’s students: 
Prof. Angia SASSI and Prof. Franco ALGOSTINO, who recalled the Master 
and Friend and the serene atmosphere of collaboration that he had been able 
to create both in teaching and in research and Eng. Dante BIGI who was the 
first Ph.D to qualify under the guidance of Prof. Cicala. 
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The Mayor of Turin, Prof. Valentino CASTELLANI and Prof. Dionigi 
GALLETTO of the University also sent their support. 

At the end of the first session, dedicated to the recollections of the 
various Organisations, colleagues and students, the meeting dedicated to the 
scientific memorials in the programme, the contents of which are given on 
the following pages, took place. 


Biografia 


Placido CICALA è nato a Messina il 9 Giugno 1910: ha frequentato il 
Politecnico di Torino, dove si è laureato in Ingegneria Meccanica nel Luglio 
1932 conseguendo la lode, ed in Ingegneria Aeronautica nel Luglio 1933, 
ancora con lode. 

Nel 1934, entrato nella carriera universitaria per suggerimento del Prof. 
M. PANETTI, è assistente ordinario; nel 1936 consegue la libera docenza in 
Meccanica Applicata, mentre già nel 1935 gli viene conferito l’incarico di 
insegnamento del corso di Calcolo di Aeromobili. Nel 1942, a soli 32 anni, è 
vincitore del concorso alla cattedra di Costruzioni Aeronautiche ed è 
nominato professore straordinario presso la facoltà di Ingegneria del 
Politecnico: nel 1945 ottiene la promozione a professore ordinario. Nel 1957 
è trasferito alla cattedra di Scienza delle Costruzioni, sempre al Politecnico 
di Torino. 

Nel periodo tra il 1948 ed il 1952, per missione del Ministero degli Esteri 
concordata con il governo argentino, il Prof. Cicala ha tenuto Corsi alla 
Escuela Superior de Aerotécnica e alla Università di Còrdoba, nonché 
presso l’Università di La Plata. 

Nel 1957 è stato visiting professor presso la Purdue University negli Stati 
Uniti, dove per 5 mesi ha partecipato a ricerche su traiettorie ottime e tenuto 
Conferenze di Calcolo delle Variazioni. 

Il Prof. Cicala ha successivamente svolto un ciclo di conferenze sulle 
strutture a guscio presso varie Università statunitensi, tra cui la Stanford 
University, la Illinois University e la Yale University. 

Fin dal 1940 è socio corrispondente della Accademia delle Scienze di 
Torino, di cui diviene socio nazionale nel 1953. Nel 1959 è socio 
corrispondente della Accademia Nazionale dei Lincei e ne diviene socio 
nazionale nel 1972. È altresì socio corrispondente dell’Istituto Lombardo, 
Accademia di Scienze e Lettere. 

Nei primi anni della carriera universitaria vinse, attraverso concorso 
nazionale per esami, un cospicuo premio, l’alunnato Gori-Feroni, 
rinnovatogli per tre anni. L’ Accademia delle Scienze di Torino gli conferì il 
premio Bressa internazionale. Fu inoltre insignito della medaglia d’oro di 
benemerenza da parte del Ministero della Pubblica Istruzione ed ottenne vari 
altri riconoscimenti. Tra questi si ricorda il conferimento del Premio Torino, 
istituito dalla Società degli Ingegneri e degli Architetti in Torino per le 
persone ed Enti che abbiano onorato con la loro attività il Piemonte. 
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Nella sessione del 1968 il premio Torino per la Classe A (opere del 
pensiero, quali contributi di studio, di ricerca e di progetto) fu conferito a 
Placido Cicala con la seguente motivazione: “Mentalità matematica degna 
delle migliori tradizioni dell'Ateneo Torinese, ha dedicato integralmente la 
sua profonda attività di studioso e di ricercatore alla soluzione dei complessi 
problemi di Aerodinamica e di Scienza delle Costruzioni nel campo della 
plasticità e sul comportamento delle strutture a parete sottile, pervenendo a 
risultati che gli conferiscono chiara fama in campo internazionale. 

In particolare si deve al Cicala la prima soluzione completa delle forze 
sull’ala oscillante ad allungamento finito che destò particolare interesse 
negli ambienti scientifici specializzati come Funzione del Cicala. 

Per la soluzione dei problemi di ottimizzazione nella meccanica del volo 
il Cicala estese per primo i metodi del calcolo tradizionale delle variazioni 
alle soluzioni discontinue. Nelle strutture a pareti sottili ha dato soluzioni 
nuove nei casi di instabilità elastica in campo non lineare. 

Docente di Scienza delle Costruzioni e di Costruzioni Aeronautiche nel 
Politecnico di Torino, ha svolto con successo il non facile compito di 
ambasciatore della cultura italiana all’estero, assumendosi l’onere di corsi di 
insegnamento presso varie università americane, dove è deferentemente 
apprezzato”. 

Il Prof. Cicala è stato nominato emerito il 13 Aprile 1987 ed ha 
continuato fino all’ultimo giorno ad occuparsi delle sue ricerche, 
completando tra l’altro la rielaborazione del Volume An Engineering 
Approach to the Calculus of Variations che è stato diffuso a cura del 
Dipartimento di Ingegneria Strutturale subito dopo la sua scomparsa, 
avvenuta a Torino il 16 Giugno 1996. 


Biography 


Placido CICALA was born in Messina on June 9" 1910. He attended the 
Politecnico di Torino where he graduated with honours in Mechanical 
Engineering in July 1932 and then in Aeronautical Engineering in July 1933, 
again with honours. 

In 1934, on the suggestion of Prof. M. PANETTI, he started his university 
career as an assistant; in 1936 he obtained the university teaching 
qualification in Meccanica Applicata, while he had been given the 
responsibility of teaching the course in Calcolo di Aeromobili already back 
in 1935. In 1942, only 32 years old, he won the chair in Costruzioni 
Aeronautiche and was named professor at the faculty of Engineering at the 
Politecnico. He was named full professor in 1945. He transferred to the 
chair of Scienza delle Costruzioni, still at the Politecnico di Torino, in 1957. 

In the period between 1948 and 1952, as an envoy of the Foreign 
Ministry and in agreement with the Argentinean Government, Prof. Cicala 
held courses at the Escuela Superior de Aerotécnica, at the University of 
Còrdoba and at the University of La Plata. 

As visiting professor at Purdue University in The United States in 1957, 
he worked for 5 months on research in optimal trajectories and held 
conferences on Variation Calculus. 

Prof. Cicala then held a series of conferences on monocoques at various 
American universities, including Stanford, Illinois and Yale. 

He became a correspondent member of the Accademia delle Scienze di 
Torino in 1940 and a national member in 1953. He became a correspondent 
member of the Accademia Nazionale dei Lincei in 1959 and a national 
member in 1972. He was also a correspondent member of the Istituto 
Lombardo, Accademia di Scienze e Lettere. 

During the early years of his university career he won the distinguished 
Gori-Feroni apprenticeship prize, for three years. The Accademia delle 
Scienze di Torino awarded him the international Bressa prize. He was also 
decorated with the Gold Medal of Merit by the Ministry of Public 
Instruction and obtained various other acknowledgements. Among many 
mention can be made of the Premio Torino, established by the Società degli 
Ingegneri e degli Architetti in Turin, for those people and organisations that 
have, through their work, brought honour to Piedmont. 

The Premio Torino-Class A  (intellectual achievements such as 
contributions through study, research and planning) was conferred on 
Placido Cicala for the following reasons: with a “mathematical mentality 
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worthy of the finest traditions of Turin University, he completely dedicated 
his profound activity as a scholar and researcher to the solving of the 
complex problems connected to Aerodynamics and Structural Engineering 
of the Construction in the field of plasticity and on the behaviour of thin 
walled structures attaining results that bestowed him with great fame in the 
international field. 

In particular mention can be made of the first complete solution of forces 
on vibrating wings of finite elongation, known as the Cicala Function, which 
created great interest in specialised scientific circles. 

Prof. Cicala was the first to extend traditional calculation methods for the 
variations of discontinuous solutions to the problems connected to the 
optimisation of flight mechanics. In thin wall structures, he produced new 
solutions in the cases of elastic instability in the non linear field. 

As a lecturer in Costruzioni Aeronautiche and Scienza delle Costruzioni 
at the Politecnico di Torino, he carried out the not easy Job of ambassador of 
Italian culture abroad with great success, taking on the burden of teaching 
courses at various American universities where he was greatly appreciated. 

Prof. Cicala was named emeritus professor on 13° April 1987 and he 
continued, up to the last, to be involved in his research programmes and 
completed, amongst others, the re-drawing up of the volume entitled An 
Engineering Approach to the Calculus of Variations which was distributed 
by the Dipartimento di Ingegneria Strutturale immediately after his death 
which occurred on June 16" 1996. 
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List of publications of Placido Cicala 


Legenda 

A. = Atti dell’Accademia delle Scienze di Torino. Serie Fisica 
Ae. = L’Aeronautica 

L. = Rendiconti dell’ Accademia Nazionale dei Lincei 


R. = Revista de la Facultad de ciencias exactas fisicas y naturales, Universidad de Cordoba 





[1] La stabilità nelle regolazioni di velocità e pressione nelle turbine a vapore, A.,Vol. 68, March 
1933, 326-336. 


[2 


La regolazione della turbina a vapore a ricupero parziale, Ricerche d’Ingegneria, year II, n. 1, 
January 1934, 22-28. 


[3 


Comportamento delle ogive coniche a velocità ultrasonore, La Ricerca Scientifica, year V, vol. II, 
n. 3-4, August 1934, 53-61. 


[4 


La flessione delle travi con piastra sottile, A., vol. 69, January 1934, 171-187, 
The bending of beams with thin tension flanges, N.A.C.A. Tec. Mem., No. 769. 


[5 


Sulla torsione delle ali a sbalzo, A., vol. 69, March 1934, pp. 313-326, 
Sulla torsione delle ali a sbalzo, Ae, vol. XIV, n. 7, July 1934, pp. 827-853. 


[6 


Schema di calcolo di un castello motore, Ae., vol. XV, n. 2, February 1935, 170-179. 
[7] La torsione dei solidi cilindrici a sezione allungata, A.,vol. 70, March 1935, 337-355. 


[8 


Il centro di taglio nei solidi cilindrici, A., vol. 70, March 1935, 356-371. 


[9 


Le azioni aerodinamiche sui profili di ala oscillanti in presenza di corrente uniforme, Memorie 
R. Acc. delle Scienze di Torino, t. 68, ser. 2, parte I, July 1935, 75-98. 


[10] Misure d'elasticità coll’uso dei pendoli accoppiati, Ricerche d’Ingegneria, year IV, n. 1, January 
1936, 1-8. 


[11] Ricerche sperimentali sulle vibrazioni flesso-torsionali di un modello di ala con rigidezze 
variabili, A., vol. 71, June 1936, 554-579. 


[12] Le azioni aerodinamiche sul profilo oscillante; Ae., vol. XVI, n.8, August 1936, 652-665. 


[13] Le oscillazioni flessotorsionali di un'ala in corrente uniforme, Ae., vol. XVI, n. 11, November 
1936, 785-801. 
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[14] Le tensioni normali nella torsione di strutture con diaframmi rigidi, Ae., vol. XVII, n.2, February 
1937, 80-119. 


[15] Ricerche sperimentali sulle azioni aerodinamiche sopra l'ala oscillante, Ae., vol. XVII, n. 5, 
May 1937, 405-414. 


[16] Sul moto non stazionario di un'ala di allungamento finito, L., vol. XXV, serie VI, August 1937, 
97-102. 


[17] Methoden zur Messung der Luftkréfte auf schwingenden Fliigeln, Vortrige der Haupt- 
versammlung d. Lilienthal Ges., 1937, 256-257. 


[18] 2 problema aerodinamico del volo ad ala battente, Ae., vol. XVII, n. 11, November 1937, 955- 
960. 


[19] Ricerche sperimentali sulle azioni aerodinamiche sopra l'ala oscillante, serie II, Ae., vol. XVII, 
n. 12, December 1937, 1043-1046. 


[20] La teoria e l'esperimento nel fenomeno delle vibrazioni alari, Ae., vol. XVIII, n. 4, April 1938, 
pp. 412-434. 
Comparison of theory with experiment in the phenomenon of wing flutter, N.A.C.A. Tec. Mem., 


n. 887. 


[21] Su/ calcolo dell'ala bilongherone con rivestimento resistente al taglio, Ae., vol. XIX, n. 1, 
January 1939, pp. 3-26. 


[22] Le oscillazioni proprie di un corpo rigido sostenuto elasticamente, A., vol. 74, April 1939, 381- 
391. 


[23] Sulle travi di altezza variabile, A., vol. 74, May 1939, 392-402. 


[24] Sul calcolo dell'ala bilongherone con rivestimento resistente al taglio, “Atti del V Convegno di 
Aerotecnica”, Napoli 1938. 


[25] Sulla stabilità dell'equilibrio elastico, A. vol. 75, January 1940, 185-222. 

[26] La stabilità elastica del cilindro in parete sottile, Ae. vol. XX, n.5, May 1940, 355-373. 

[27] Sul calcolo dei solidi cilindrici in parete sottile. Parte I, Ae., vol. XX, n. 6, June 1940, 460-479. 
[28] Su/ calcolo dei solidi cilindrici in parete sottile. Parte II, Ae., vol. XX, n. 7, July 1940, 595-612. 


[29] Su! calcolo dei solidi cilindrici in parete sottile. Parte III, Ae., vol. XX, n. 10, October 1940, 
735-746. 


[30] Sul calcolo delle strutture alari a guscio, “Atti del VI Convegno di Aerotecnica”, June 1940. 


[31] Sulla teoria non lineare di elasticità, A., vol. 76, December 1940, 94-104. 
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[32] Sui cassoni alari con pareti a traliccio diagonale, Ae., vol. XXI, n. 3, March 1941, 183-204. 


[33] Le teorie approssimate dell'ala oscillante di allungamento finito, A. vol. 76, March 1941, 389- 
416. 


[34] Ricerche sperimentali sulle azioni aerodinamiche sopra l'ala oscillante, serie III, Ae., vol. XXI, 
n. 1, January 1941, 46-53. 


[35] Lo stato attuale delle ricerche sul moto instazionario di una superficie portante. Parte I - Il 
problema bidimensionale, Ae., vol. XXI, n. 9, September 1941, 557-591. 


[36] Idem. Parte Il - L'ala di allungamento finito, Ae., vol. XXI, n. 11, November 1941, 671-685. 

[37] Idem. Parte III - Resistenza e propulsione, Ae., vol. XXI, n. 12, December 1941, 759-769. 

[38] Idem. Parte IV - Le ricerche sperimentali, Ae., vol. XXI, n. 12, December 1941, 770-773. 
L'état actuel des recherches sur le mouvement non-stationnaire d'une surface portante, 
O.N.E.R.A., Trad. Tech., n. 1067. 

[39] Sul calcolo delle strutture a cassone, A., vol. 77, March 1942, 209-241. 

[40] Sul calcolo del guscio a quattro correnti, Ae., vol. XXII, n. 12, December 1942, pp. 588-621. 


[40bis | Calcolo delle Strutture di Aeromobili, Edit. V. Giorgio, Torino, 1943, 532 p. 


[41] 1! cilindro in parete sottile compresso assialmente. Nuovo orientamento dell'indagine sulla 
stabilità elastica, Ae., vol. XXIV, June 1944, 3-18. 


[42] La struttura a guscio soggetta a carico trasversale uniforme, A., vol. 80, June 1945, 114-120. 


[43] Sul dimensionamento dei correnti principali nelle strutture a guscio, A., vol. 80, June 1945, 121- 
109; 


[44] L'interruzione di parete nelle strutture alari a guscio, A., vol. 80, June 1945, 128-133. 
[45] Effetti di un carico tagliante nelle strutture alari a guscio, A., vol. 80, June 1945, 134-139. 
[46] Le tensioni caratteristiche per il guscio cilindrico, L., serie VIII, vol. I, n. 6, June 1946, 735-739. 


[47] Lo stato di tensione nel guscio cilindrico presso una accidentalità, L., ser. VII, vol. I, n. 10, 
October 1946, 1069-1073. 


[48] Sul calcolo delle strutture a guscio, Ae., Vol. XXVI, n. 3, September 1946, 138-148. Ae., Vol. 
XXVI, n. 4, December 1946, 216-229. Ae., Vol. XXVII, n.1, February 1947, 31-43. Ae., Vol. 
XXVII, n. 2, April 1947, 133-146. 


[49] Recenti ricerche sulle vibrazioni alari, Ae., vol. XXVII, n. 1, February 1947, 54-59. 
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[50] Recenti ricerche sulle strutture a guscio, Ae., vol. XXVII, n. 3, June 1947, 230-235. 


[51] /! guscio a direttrice circolare, con correnti e pannelli uguali, soggetto a carichi trasversali, 
Monogr. Scientif. Ministero Areonautica, n. 7, August 1947. 


[52] Sull'analisi delle piccole deformazioni in campo elasto-plastico, L., ser. VII, vol. II, July 1947, 
325-329. 


[53] Effects of cutouts in semimonocoque structures, Journ. Aeron. Sciences, vol. 15, n. 3, March 
1948, 171-179. 


[54] Problemas de calculo de las estructuras en pared delgada, R, aîio XII, n. 1-2, 1949, 165-193. 


[55] Efectos de cargas torsoras actuando cerca del fuselaje sobre estructuras de alas con dos 
largueros, R., afio XII, n. 3, 1949, 517-554. 


[56] Sobre las vibraciones alares de flexiòn y torsiòn, R., aio XII, n.4, 1949, 925-951. 


[57] Sul comportamento elastico di una parete sottile quasi cilndrica, L., ser. VII, vol. VII, n. 1-4, 
July 1949, 99-103. 


[58] Sulle deformazioni plastiche, L., ser. VII, vol. VIII, n. 8, May 1950, 583-586. 


[59] On plastic buckling of a compressed strip, Journ. Aeron. Sciences, vol. 17, n. 6, June 1950, 378- 
379. 


[60] Su/ carico critico di una piastra compressa oltre il limite elastico, L., ser. VII, vol. IX, n. 1-2, 
July 1950, pp. 67-71. 


[61] Sobre la teorià de Batdorf y Budianski de la deformaciòn plastica, R., aîio XIII, n.2, 1950, 401- 
415. 


[62] Column buckling in the elastoplastic range, Journ. Aeron. Sciences, vol. 17, n. 8, August 1950, 
508-512. 


[63] La distribuzione ottima degli spessori del rivestimento alare, A., vol. 85, November 1950, 143- 
158. 


[64] Sul criterio di rigidezza delle strutture alari, Termotecnica, vol. 5, n. 2, February 1951, 45-49, 


[65] Sobre un problema de pandeo de placas en campo elasto-plastico, R., aîio XIV, n. 1, 1951, 29- 
45. 


[66] Procedimiento aproximado para el estudio del pandeo elastoplastico, R., aîîo XIV, n. 1, 1951, 
62-78. 


[67] On plastic buckling of plates and a theory of plastic slip, Journ. Aeron. Science, vol. 18, n. 4, 
April 1951, 285-286. 
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[68] The effect of initial deformations on the behaviour of a cylindrical shell under axial 
compression, Quarterly of Applied Mathematics, vol. IX, n. 3, October 1951, pp. 273-293. 


[69] Sulla determinazione del comportamento del sistema pneumatico ammortizzatore, A., vol. 86, 
December 1951, 25-35. 


[70] Determination of frequencies and modes above the fundamental by matrix iteration, Journ. 
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Modernity and potential developments 
of Placido Cicala’s thought 


Extended text of the opening address to the 
Memorial Day of Placido Cicala 


Ettore ANTONA 


1. Introduction 


“The greatness of a scientific contribution does not mainly lie in the 
exactness of the affirmations but rather in its capability to induce 
developments and innovative ideas.” Such affirmation, considered valid by 
many eminent scientists, can be undertaken as a guide to this text, which 
would demonstrate some of the potentialities of the work of Placido Cicala. 
His work is enormous and includes diversified and numerous contributions 
in various directions that are inter-connected by the mutual origin of an 
aeronautical-structural sector of mathematical physics and, above all, by a 
vision that characterised Cicala, about the ways through which one could 
contribute to the progress in knowledge. To be precise, it should also be 
mentioned that such “exactness”, which, according to the mentioned saying, 
is not strictly necessary for the greatness of a scientific contribution, is very 
much present in the whole work of Cicala at a high degree, for there are no 
formal errors, no incoherence with respect to the assumptions and, finally, 
there is a presence of a great deal of information on other scientists’ 
contributions. Cicala’s writings from this point of view, are examples 
because of their exactitude and correctness. This is obviously a merit, 
however we will devote ourselves to examining the potentiality of his 
thought. To those who would object that, at the moment, Cicala’s idea do 
not seem to have a large following, it is possible to answer that: 

- ideas must be known in order to be developed, and further efforts are 
necessary in order to present Cicala’s ideas to a wider number of 
scientists; 

- for many aspects the mass of the potential users of Cicala’s ideas might 
need more time; 
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- the rapid transformation in progress relating to applicative analyses, 
caused by the spreading of digital elaborators and their connected 
philosophies, makes more difficult the understanding of the potentiality 
of Cicala’s ideas and their applicability to the new realities. 

Let°s now consider the guidelines of the thinking on modernity and on the 
potential developments of Cicala’s ideas. 

He didn’t like the too “generalistic” formulation of problems, I do not know 
by means of which combination of instinet and rationality preferring to face 
well defined questions, whose treatment makes it possible to examine the 
capability of a theory to come close to the exactness. He often said: “I°m not 
interested in generalisations; I prefer to face a well defined problem.” This 
did not stop him, however, from inserting that problem into a context that 
could be understood from all points of view and for which the solution was 
then found. 

When facing a problem he was always well aware of the importance of 
inserting a solution or an analysis into a growing approximation treatment 
context, with the possibility of estimating the obtained approximations. 

This coherent and convincing outlook had a great influence on me, so 
much that, in a more general environment, I found myself thinking, without 
obviously meaning to put on others the responsability of my mistakes on 
that route, that one can find starting points for further steps in the 
advancement of knowledge, in a re-evaluation of the way of operating of 
technical sciences (of which Cicala was one of the greatest spokesmen) 
which accepts approximations and the limitation of both instruments and 
analyses. 

Knowledge is not in fact usually obtained in complete form but rather in 
an approximate way. This does not mean that in some simple cases it cannot 
be in a complete form but the assumption a priori that a complete form could 
be always possible is a leading astray and can cause method errors. The 
problem is instead to establish whether the obtained approximation is 
sufficient for the end purposes which are searched for each time. 

In the same way, language in general gives approximate descriptions and 
their indefinite improvement possibility must suffice. 

Thought, even though conceiving complete concepts, should work on 
approximations that language can furnish and on the knowledge that the 
obtained approximation makes available. 

This is in contrast with the positions that maintain that those concepts 
that are dealt with by science and in particular by technical sciences are 
“pseudo concepts” that are not worthy of philosophy. Such “pseudo 
concepts” are on the contrary the only coherent ones with an acceptable 
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outlook of the knowledge, while one should not trust any concept for which 
it is impossible to identify an approximate determination. 

Engineers make great use of mathematical models, here indicated as 
“engineering models.” 

Among the other interpretations, mathematics can be considered as an 
instrument for the solution of practical problems. Mathematics in fact 
includes numerous disciplines. Many models, indicated as “mathematical 
models”, founded on axioms and rules, are necessary to formally make up 
the various disciplines. 

Mathematics also makes instruments available for the solution of 
problems deriving from humane activity and this by means of mathematical 
models that contain approximations of the effective behaviour of practical 
systems of various degrees of complication. This would seem to be a 
completely different meaning of the expression “mathematical models.” 

In all cases such models are founded on rules and axioms, including 
those that are necessary for the founding of the involved mathematical 
disciplines and those that are necessary to create a correspondence between 
the mathematical description and the physical phenomena. In other words, 
both the mathematical models and the engineering models have an axiomatic 
nature in common that would justify the unified use of the expression 
“mathematical models.” 

Nature and structure of mathematical models offer two important 
subjects for discussion. One is the problem of coherence and completeness 
of the various axioms. The other is that of the possibilty of each model to be 
included in a gradually improving process of approximation towards the 
exact solution. 

As far as formal systems for mathematics are concerned, coherence and 
completeness of the axioms were a priori presumed for a long time. There 
was however a great crisis after the fundamental contribution by Gòdel, who 
demonstrated, in 1931, that each axiomatic structure, if pretended to be 
coherent, it results to be incomplete and if pretended to be complete, it 
results to be incoherent. His demonstration made use of mathematical 
evaluations of metamathematical affirmations. Nevertheless, when 
necessary, mathematicians continue to use their formal systems, probably 
supported by the knowledge that the problem shown by Gédel is very 
remote, in terms of the number of operations, besides the deductions or 
demonstrations they require. 


"In [15], where Cicala’s work was also taken into account, the A. had already discussed the thesis 


proposed up to the end of this point making large use of such quoted references. 
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In practical activities the role of mathematicians is mainly that of 
formulating mathematical models of practical problems and not that of 
developing theories or demonstrations on the convergence of numerical 
processes. The importance of this aspect of mathematics has considerably 
increased over the last twenty years, particularly in the research and 
industrial development activities where the majority of work is done by 
engineers. 

Equations describing real problems are usually simplified with extra 
axioms or, according to Cicala, with asymptotic approaches. 

In order to resolve engineering problems, many types of functions are 
introduced: physical principles or laws, behavioural equations (simplified or 
non simplified), and simplifying hypotheses. If such functions are thought of 
as being axioms, each mathematical model that has been created for a 
practical application can be considered as a formal structure to which a 
deductive analysis can be applied. 

The axioms that introduce approximations are usually reasons for 
incoherence, in particular’ when they represent limitations or 
characterisations of the unknown. One way of avoiding such incoherence 
risks is that of reducing the set of problems covered by an axiomatic system 
to a subset of the problems for which the axioms are acceptable. 

Engineering is very much founded on approximations, in order to 
simplify equations. Approximations are introduced while less important 
physical effects are neglected and other simplifications are introduced 
through axioms concerning the solutions or other behavioural parameters. 

As far as mathematical models applied in engineering are concerned, the 
various approaches can be classified in the following manner, in particular 
after the fundamental contributions made by Cicala: I.axiomatic approaches 
that require axioms concerning the unknown functions of the problem and 
which directly resolve the mathematically obtained models; 2.aprioristic 
approaches that obtain the mathematical models through axioms concerning 
the unknown functions, but which resolve such models through asymptotic 
concepts; asymptotic approaches that obtain mathematical models and 
resolve them through asymptotic concepts, with various advantages such as 
the determination of a sequence of gradually improving approximation steps 
and the asymptotic coherence of each step. 

AII these approaches should obviously be considered heuristic. 

Another fundamental aspect is that of the mutual independence of the 
axioms, or rather of the minimum number of axioms at the basis of a 
mathematical model. This is particularly important when “principles” are 
numbered in the role of the axioms. The position of Cicala that emerges in 
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particular from his contribution to variational calculus is that of lowering all 
the properties of a minimum axiom set, at the same time formulating the 
mathematical model, in such a way as to be able to consider the peculiarities 
of the practical problems, therein included, for example, the continuous 
solutions and the piecewise solutions and the conditions posed under 
inequality form. 


2. Synthesis of Cicala’s activities 


In this moment that we are going through, there is the risk that a 
monotonous image of the progress being made in technical sciences, which 
favours automatic calculation, does not fully acknowledge the teaching that 
can be derived from Cicala’s work, even though his ideas are not in fact 
offered as an alternative to automatic calculation; on the contrary they are 
susceptible of developing their capabilities in a great deal apart from 
contributing to the fundamental theoretical clarification of the physical- 
mathematical aspects of the problems he dealt with. 

Cicala’s activities, if one wants to schematise them, were basically 
carried out in the following fields: reinforced shells, snap buckling 
instability; variational calculus; asymptotic approaches: elastoplasticity 
theory and also non homogeneous materials at elastoplastic behaviour. 

The problem of aerodynamic actions on profiles and wings in non 
stationary movement were among the first, as far as time is concerned, to be 
faced by Cicala. His research spread internationally far and wide: the 
integrodifferential equation that determines the distribution of the total 
circulation along the wing, in the case of non stationary movement, and 
which is named after him, has been compared because of its importance 
with that of Prandtl, which is also a particular case. He was a forerunner in 
this field. Many other authors came close to him and their methods, later on, 
compared by Cicala himself to his own, result to be basically equivalent. 

As far as the snap buckling phenomenon is concerned, in a fundamental 
work of 1941, Kàrman and Tsien were together ascribed the merit of having 
“rediscovered” the importance of non linear terms in the analyses of 
structures containing instabilities of such a type. The scientific community 
was then about ten years ahead of an unexplainable inadequacy of the 
theories, all because of a strange blindness founded on linearized 
approaches. Back in 1940 Placido Cicala also carefully examined linearized 
theories, enunciating their insufficiency. Later on and independently (this 
being during the war period with the consequent communication difficulties) 
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he faced the problem with non linear theories, explaining how initial 
imperfections of form could induce the experimental behaviour that shows 
critical levels of an order of a fifth of those that are justifiable with linear 
theories. Under different conditions it would have been possible to attribute 
him with a share of the priority. In following works the results obtained by 
Cicala became enriched by discoveries of the physical behaviour expected 
by the theory and, as far as I know, not used. although it is really possible, 
even on a practical level. I am referring to the observations, that are 
perfectly in tune with the completeness of horizons of his intuition, which, 
as undesired and casual geometric imperfections drastically lower the limit 
of stability, well aimed initial forms (that can be in contrast with phenomena 
connected to instability) can in fact raise the critical limits to the “perfect” 
forms, for example in the cylinder in thin walls. 

In the last period he also dedicated himself to other subjects such as non 
homogeneous materials and elastoplastic behaviour. In works that were 
outstanding because of their depth and extensiveness, he examined the 
foundations of the various known approaches with clarification and 
generalisation that, yet again, highlight his extraordinary capacity to again 
give his higher level outlook to the majority of treatments and theories. 

Space requirements make it impossible to give details on the 
aforementioned fields to illustrate their modernity, in particular for the times 
in which they were dealt with, and their possibility of development. It is 
therefore necessary to limit the aim of this text to the remaining fields: 
reinforced shells, variational calculus and asymptotic approaches. On the 
asymptotic approaches, because of the amplitude of the original 
contribution, our exposition will be particularly brief, in fact it will be 
restricted to a review of Cicala’s results and of their fundamental concepts. 


3. Characteristic states 
3.1 Reinforced shells 


On the linear theory of reinforced shells, an item of structural analysis, 
the research of Placido Cicala, according to Argyris, has the same level as 
that of the best of German literature. Apart from these recognitions, Cicala, 
with original intuition and rigorous procedures, examines and goes to a great 
depth into those aspects that are not so frequent in international literature, 
such as: the local stress states due to constraints and discontinuity and the 
effects of the deformability of transversal stiffening; the existence of 
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characteristic stress and strain states, which can constitute complete 
orthogonal bases for the development in summation of the solutions, which 
was a completely original subject due to Cicala. 

It should also be noted that, even the effects were studied to some extent 
by other authors, they were inscribed by Cicala in a stimulating vision of 
treatment sequences that were capable of growing approximation. The idea, 
completely original, of the “characteristic states” in particular dealt with 
here, favours, above all, transpositions to other approaches such as finite 
elements, with the advantage of clarity and power of analysis. 

Cicala elaborated the idea of the “characteristic states” in the case of 
“cylindrical” reinforced shells with only cuts in panels and made up of 
identical modules (delimited by transversal stiffness) all with the same 
length. 

The configuration of stresses and displacements that is determined on the 
basis of the elementary theory respects the continuity of the body and the 
equilibrium conditions on the whole structure, with the exception of those 
sections where there exist casualness, represented by external loads, 
constraints or variations in the dimensions of a structural element. 

It is possible to define corrective states as the states that should be added 
to those of the elementary theory to compensate the discontinuity of the 
stresses and eliminate the congruence defects. 

Corrective stresses, whether normal or tangential stresses, are 
conventionally marked with an asterisk in the following pages. 

As the stresses of the elementary theory balance the applied loads, in 
each section, the corrective stresses should therefore have a zero resultant 
and a zero resultant moment, in order not to vary the already present global 
(but not local) equilibrium state. 

As far as linear approximations are concerned it is possible to 
demonstrate that corrective stresses in shells of constant dimensions, can be 
expressed as linear combinations of characteristic stress systems. 

A system of characteristic stresses (a system of stresses having 
characteristic propagation) is a system of stresses, which, applied to 
whatsoever section of a shell, composed of ideals modules, induces 
proportional stresses in the other sections of the shell, according to a 
geometric series. 


er] 


x 
“ Panels are longitudinally stiffened by stringers and are transversally stiffened by ribs (wings, tails) or 
frames (fuselage, nacelles). Here the word “frame” is used for any transverse stiffener. 
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Figura 3.1 
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Figura 3.2 


The situation that is found in the various sections between modules (the 
initial section has the order number equal to zero) is the following 


[S]: Si c AS, 0 = Z% ’ S;, = 45,1 = À Sj ’ S;3 45, i ÀS; ’ 


where S means a longitudinal stress in a stringer. As congruence between 
mudule and module is respected, the displacements of the separation 
sections obviously result to be proportional in geometrie series, [5]; also 
shear stresses result to be proportional in geometrie series: 


li = Ao ro Mi, I, = Mi, "I Ati ? l,3 Gi A, CAL 


I 
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The 4 factor is the relationship of the geometrie series. In other terms [5], 


(3.1) 
=A;,. (3.2) 


An example of characteristic propagation is given (with Z</) in figures 
3.1 and 3.2. 
The number of characteristic stress systems for the closed shell (twice 


connected), with m currents, is p = 2 (m-3). 

This number results to be halved if deformabilities of the frames are 
neglected. If the section is m’#2 times connected, presenting nm’ longitudinal 
frames, m/ characteristic systems are added that disappear in the case of 
rigid transverse stiffness. 

From such equations only the stress systems that have a zero resultants 
and zero resultant moments, called autoequilibrated, are searched for, 
otherwise the equilibrium of the reinforced shell, would not be possible. The 
value of Z=/ in such a case, occurs for rigid movements of the structure in 
which the stresses are null. By not imposing autoequilibrium one could also 
obtains the solutions relative to the states that propagate (statically) without 
changes at every distance. The corrective stresses can be expressed through 
a linear combination of characteristic propagation systems, [5]. 


P 
Top Mo 07 PORT (3.3) 
j=l 


a, = coefficients of the linear combination (a, depend on the form and 


intensity of the singularities that generate the corrective stresses), 
A = the characteristic stress of the j-th system, present on the i-th stringer, 


ine = the coefficient of the j-th system, which expresses the attenuation 


that such a system undergoes in the r-th section, in comparison to the section 
in which the singularities can be found, 


07 ,= the corrective stress of the i-th stringer of the r-th section. 


Eq. (3.3) presents the corrective stresses as perturbations originating in 
the section where the inducing singularity is present; therefore, for the 
calculation of the corrective stresses in each section it is necessary to take 
into consideration all the singularities present in the structure. 
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In practise, for the attenuation of the characteristic stresses it is possible 
to take into consideration, without great errors, only the casualness that are 
closest to the section under examination. 

When one considers that the system of frames is continuous (see also 


section 3.2), that is, that the ordinates are infinitely close, the Ars factor 


present in the summation is substituted by an exponential of the type e !”, 
where p, is a typical coefficient of the characteristic system, and 2 the 
distance from the singularities, referring to a reference length. 

The research of the characteristic systems is carried out by taking into 
account the reinforced shell laws and the property of propagation of the 
characteristic systems. 

A segment of a structure of constant dimensions is considered which is 
defined by two frames, between which other frames are interposed that are 
equal and equidistant. 

With a larger analysis than that proposed by Cicala, who considered the 
case in which on each transversal section the applied forces are in 
equilibrium (self-equilibrated), each box must however result to be in 
equilibrium. The following equations, to which type S forces do not 
contribute, are therefore valid 


Yi, -x,)1-2)=0. (3.4) 
Y1,0, TS; )1-4)=0, (3.5) 
Y,,9(1-2) =0. (3.6) 


Such equations are verified if Z =1, that is, if the system of forces applied to 
the extreme of the module cross it unaltered. In the opposite case (4 # 1) the 
following equations occur: 


n 3.4' 
Di (a, da X;-1 ) = 0, i 
i=l 
a (3.539) 
deboli, hi Vi ) 77 0, 
ial 

Std at 


The remaining equilibrium equations of the module: 
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) S..(1-4)=0, 
Lav) 3.7) 
LP F Do — XK) (3.8) 
S..y;(1-4)= VExi;ee 
L AI ) Zt50 Via (3.9) 


are not independent from (3.4), (3.5) and (3.6). If Z =I, (3.7) is 
automatically verified but (3.8) and (3.9) are not verified and it is also 
necessary to verify (3.4’) and (3.5°), that is, the zeroing of the components 
in x and y of the shear on each section. The only states that can pass a box 
unaltered (4 =1) are therefore constant longitudinal forces and bending 
moments. 

In the case of 7 #1, (3.4°), (3.5’) and (3.6’), considered by Cicala, are 
valid and (3.7). (3.8) and (3.9) reduce to three equations that can be obtained 
by the three previous ones and generate selfequilibrated systems on each 
section. 

For characteristic systems that are equilibrated on each section (4 #1), 
going back more closely to Cicala’s treatment, the following relations exist 
[5]: 

S 


t 


=7S,, (3.10) 


i,r+l i,r 


=Af,, (3.11) 


hral Cla 

where 4 is the propagation parameter having a value that is constant on 
the whole structure, whatever could be the node (i, ) or the panel (i, r) . 

For a system of characteristic forces it is therefore possible to write, by 
using the last relations, [5]: 





1-4 
ga DI 3912 
GIA i,r i+1,r i,r ( ) 
and then, [5]: 
mic; j VERSI TVirti -2v,,)= 
I 
4+4+— 
; Sora É 
pprranvente itadeniotitta wé Biz grgpi datephinrAbopre (3.13) 
6E A, A; PARI E; i 


The frames are loaded by the differences between the flows in the same 
panel on one hand and by the ordinates on the other, that is, by the 
quantities: 
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La db » (4 — I): li, r (3. 14) 
On the basis of all this one can write, [5]: 


CV Rig — (x, — i Xi cata -(y, - vt, 20,0, 


= 7 1). [0'0,dv 
where o; is the stress caused by the tangential distribution with the density 
t, between the i-th and its preceding currents and Where the integrals are 
extended to the whole volume of the z-th frame. 

The corresponding expressions for the ( r-/) and ( r+/) frames have the 


(3.15) 


I 
second member multiplied respectively by Pa and by 7. 


By writing these two expressions and joining them with that written for 
the ( r) frame, the following relations are obtained [5]: 


Cc; 1 (V,.., ni Vizi LR 2v), )- (x, (di Xj-1 ). (60% Py 26) Gio Ou )+ 
Tr (y, To Vi ). (22, wi 2, ii 7,41 )- 2A, (9... mi 20 + Bua )= (3, 16) 


-(# =3AR3 -1) [0-'0,dv. 
À Q 


The following quantities: 


VV ep ZVpprkiMoas 


< p> as xe 26, x Sri 9 

M =M72,+M > 

0 = 0,; a 20, + 9,,ì 
are proportional, as already mentioned, to the derivative, made on an 
interval with a constant step, of the size to which reference is made. 

Let's suppose that the section has m wall arches and presents m’ internal 
panels, i.e. it is 2+m' times connected. From equations of the (3.12) type it 
is possible to write as many equations as there are stringers, that is, n. The 
(3.13) equation can be written as many times as the wall span, that is, m+m?”. 
The equation (3.16) can be written the same number of times, m+m’; there 


are also the three equations (3.4,5,6). There is therefore a total of 
3m+2m'+3 available equations. 


i,r-l 


(3.17) 
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The unknowns are: the normal S forces in number of m; the shear flows 
t,, in numbers of m+m’; the quantity Vi, in number of m+m’; and the three 
quantities é'’, 777, 0”. 

The number of equations equals the number of unknowns. The stresses in 
the frames depend on the flows that act on their boundary; by supposing that 
it is possible to express the integrals according to the member of (3.16) in 
function of the /,, flows, the system of equations that is obtained is a 
homogeneous system of (3m-+2m +3) equations in (3m+2m +3) unknowns. 

The condition where the determinant of its coefficients is zero is an 
equation in Z that can be reduced to an equation in the unknown: 


pipi er 
7 


which appears in it with exponents from a maximum of 2(m-3+m') to a 
minimum of m”: 


FO agi Si 


By excluding the solution 4 = 0 (movement of rigid body), the degree of the 
polynomial can be lowered to m'+2m-6. In the case where the deformability 
of the frames is neglected, it descends to m'+m-3. 

Thus one arrives at the conclusion that there exists (referring to the case 


of deformable frames) a total of m’+2m-6 solutions 4 k of the 


characteristic. polynomial. For each solution 4 » there are two 
corresponding solutions 4, each of which is the inverse of the other. For 


each of the m'+2m-6 solutions A .. one determines, by resolving the 


h» 
homogeneous system that has reduced the number of its equations to 
m’+2m-6, the eigenvector, exception made for a constant factor, which 
constitutes a characteristic system. 

The number of characteristic systems that is thus constructed is m’'+2m- 
6. The eigenvector can be composed of m+m'’-3 shear forces and of m-3 
normal forces. The 6 static equilibrium forces allow one to calculate the 6 
normal forces. The characteristic. systems are obviously among the 
orthogonal ones and allow one to obtain, by linear combinations, any 
selfequilibrated state applied to a section. 


In the case of infinitely rigid frames, the number of A and therefore of 


the characteristic systems that one determines is m-3, and the values of 4 
result to be real and positive. 
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Otherwise, the Z can result to be complex and therefore suitable for 
complex conjugated pairs and in the same time for pairs of values and their 
inverse. The complex ones therefore result to be grouped in set of four 
values inside which the two mentioned properties are present. If m'+2m-6 is 


odd, at least one of the 4 results to be real, generating at least two real 4. It 
is interesting to note that, in the proposed formula, the problem can always 


be reduced to the determination of 4 and they may or may not have 
complex 4. 

The displacements given by a system of characteristic stresses can now 
be analysed, by first of all considering the “self-equilibrated” systems 

Let’s consider a reinforced shell structure that extends from the frame 
number zero to infinity, towards the positive z coordinate. Let’s suppose that 
such a structure is subject, in correspondence to frame number zero, to a set 
of characteristic loads, corresponding to a certain factor of propagation 4; 
#1, in the following frames number r the normal stresses are given by: 


O;, =: À;, 5 Gi n À;r i O ;.; (3.18) 
where (o 38 are the characteristic stresses of the j-th system and on the i-th 


stringer, that are applied in correspondence to the frames numbers zero. 
In the panels included between the frames number zero and one the 
flows are: 


bd (3.19) 
where Ù, is the characteristic flow of the j-th system and on the i-th 


span, applied in correspondence to the frame number zero. 
In the following panels of each span, the flow / varies according to a 


geometric progression of À ; ratio. From the two parts of the frame number r 
there will therefore be flows whose difference is: 
trogiot,g 4A, bet} (3.20) 
If the stresses inside the frame are indicated with do when the ordinate is 
loaded on its boundary by te flows, it occurs that, in the following frames, 
the o stresses have the values given by: 

Gy =A;(4,-1D:0,,;. (3.21) 


This occurs as each frame is loaded by the differences that exist between the 
flows of the same span, in the two panels separated by the frame itself. 
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Two distinct modules result to be subject to stress data at the terminals 
that are part of the 4% relation, if they are distanced from other m-l 


modules. Their deformation will be, in a linear behaviour, in the same À 


relation. In particular, each module does not undergo, under each 
characteristic system of 4, eigenvalues, either variation in length or 


rotations relative to the terminal sections. This can easily be verified with 
elementary applications of the virtual works. When referring therefore to the 
displacements of the elements of the terminal sections to the positions in the 
unloaded box it can be concluded that, for selfequilibrated systems, not only 
the stresses at the terminal sections but also the displacements are within the 


A; relation. 


With an increase in r, because the 4; module is lower than unity, all the 
stresses tend to disappear as well as do the deformations of the structural 
elements. The displacements of the far sections (z tending towards infinity) 
in the axial direction are therefore null. In such conditions, the axial 
displacement w,, of the (i ,r) point is also represented by the integral of the 


, ; î o GL 5 
section number r extended towards infinity , by the E strain in the i-th 


stringer. By taking into consideration the initial values of the stresses and 
the fact that they vary linearly between one frame and another, one obtains: 


a-(1+4,) a 
Wi} su AD ji 
AM 0 ao 


The systems connected to the Z = / value, in turn, generate  symmetric 
displacements in each module if they refer to the middle section ( with 
respect to which the structure is symmetric and the stresses are symmetric. 

When wishing to identify a condition of unitary ratio (4 =1) even 
between the displacements of the terminal sections, it is possible to think to 
infinite displacements (in translation and/or rotation) of the module, which 
are added to those that were previously considered. 

It should be noted, apart from what was stated by Cicala, that in the case 
of real 4, they can , at least in part, result to be negative. This depends on 
the approximation that goes under the name “reinforced shell”, when the 
structure is sub-divided in length into “finite” modules. In the following 
paragraphs, in which an approach for infinitesimal modules is considered, it 
is possible to see that in a finite span, each 4, if real, is positive. 


(3.22) 
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3.2 Propagation by shear 


In order to show the ductility of the idea, it is possible to use an approach 
that permits the propagation by shear of the panels to identify characteristic 
states that propagate with constant attenuation on constant segments. A 
continuous frames system is (implicitly) supposed in these approaches. 

One has, if m is the number of stringers, 2 - 6 equations in 2m - 6 
unknowns for each section between two modules. If the frames are 
deformable, they, in order to be in equilibrium under shear flow differences 
transmitted by the shell, deform on their plane. It is 
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m =7= number of ribs and skin panels 
m° =2= number of internal panels 





MARK ELEMENT NUMBER 


sigrl stringer m 


AE internal panel m' 


Miao skin panel m 


lr 4° frame 


Fig. 3.3 
therefore necessary to consider the displacements of the nodes in such a 
plane: it is sufficient to consider those acting as the directrix of each shell 
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span. This introduces, if m-+m’ is the number of the spans, other 2(#-4m/) - 6 
equations in 2(#+m’) - 6 unknowns. In many cases it is justifiable to neglect 
the deformability of the frames (see Fig.3.3). 

The autoequilibrated force systems that should be introduced at the 
terminals of the modules can be obtained as linear combinations of 
autoequilibrated characteristic systems, which have the property of 
propagating in the structure, each one attenuating according to its own law, 
in such a way to be however represented, at least by constants, by the same 
graph on each transversal section. 





3 


= 


Fig3 4 


Such characteristic systems are equal in number to that of the parameters 
that define the stresses on the transversal section, that is, 272+m 6 and every 
one of them is orthogonal to each other. 
As on each cut two series of characteristic systems are introduced, one for 
each module, for each cut there are 27+m 6 systems, as many as there are 
the unknowns of the problem concerning the corrective stresses. 

Let's consider, as an example, a transversal section of a cylindrical 
section with constant dimensions, composed of m stringers and 2m+m "6 
panels forming a closed section, with infinitely rigid and infinitely close 
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frames'. As the transversal sections do not undergo deformations in their 
plane, the state of deformation of a module is completely described by the 
displacements (see Fig. 3.4). 

The following relation are valid: 


dw dF 
F.=ES, —, Ci-V g,, 
Mr e LU 


W-wW, dé 
= Gs sitana Sagl a do, ate pa ; 
Gi, { c; LA UV Vi dz cu) 


dF, dF, 
dali Dig her 2 Z4%y Dar =) Budgzo 
pa Ì LI pi 


(3.23) 





di dq, 
pa ina) Lita > Xi x;=0, Lodi Pi V;= 0, 
IPA 


Ap;F,+bjly;= € +8, 4 +15 (p=1,2,3) (3.26) 
indicating by w;(z) the displacement of each stringer section in the 
longitudinal direction and by Az), é(2), 7(2), respectively the rotation and 
translation according to two axes of the frame in its plane. The stress state 
involves the forces F;(z) in the stringer and the flows g;;(z) in the panels. The 
(3.5) equations are obviously valid. The constants ;, 8}, and y;,, obviously 
depend on the geometry of the transversal section and on the choice of the x 
and y axes. 

As the module is only loaded at its extremes, the equilibrium equations of a 
interval dz of the module itself are valid (3.24) and (3.25). 

Three other equations are obtained by expressing the congruence 
between the deformation of the stringers and panels and the unknowns 9 é 
mn. This can be done through the virtual work equation, using three 
equilibrated conditions that are not linearly dependent between each other 
and have the moment resultant not null. Such conditions, applied to an 
infinitesimal dz, work on the effective displacements (congruent). Three 
relations of the (3.26) type are thus obtained. 


3 The A. has already discussed in [6] the argument reported up to the end of this point, making large use 
of the quoted reference. 
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The (3.23), (3.24), (3.25) and (3.26) equations make up a system of 
3m+m'"+9 homogeneous linear equations in 3m2+7m3 unknown functions. 


By searching for exponential solutions f = f, ‘e a matrix is obtained of 


the coefficients in which 7 lines must result to be linear combinations of the 
previous ones. For 6 of these, which reduce the number of equations to that 
of the unknowns, this is true independently of the value of 4, in that the 6 
equilibrium equations (3.6) and (3.7) are also implicitly contained in (3.5). It 
is therefore sufficient to select 3+m 3 among the written equations. In 
order to obtain that the determinant of the coefficients of the so formed 
system is null, an equation is placed in 4 in which, as can be verified by the 
opportune choice of available equations to form the system, the maximum 
power of 4 is 2m1+3 and the minimum m+6. The m+6 zero roots correspond 
to solutions such as translation or rigid rotation, uniform compression, 
uniform bending or they should be rejected. The not zero m-3 solutions 
correspond to autosolutions in which, as can easily be verified, the forces 
acting on each section are autoequilibrated. These also result to be 
orthogonal between themselves. 

The number of characteristic states that are identified is identical to that 
obtained with discrete frames, that is, with the real reinforced shell. 


NUMBER OF THE AUTOEQUILIBRATED CHARACTERISTIC SYSTEMS 


FRAMES DISCRETE CONTINUOUS 
' 


3.3 Characteristic states in fem approaches. 











3.3.1 General considerations 


The application of the idea of the approaches to finite elements has 
proved to be fruitful, allowing the identification of states, stresses and 
displacements on transversal sections, that on constant lengths attenuate in a 
constant way, with determination of the same states (eigenvector) and of the 
relative attenuations (eigenvalues) It has been possible to isolate the 


MODERNITY OF PLACIDO CICALA’S THOUGHT 49 


movements (with unloaded structure) of rigid bodies and show the effect of 
the transversal deformability that complex autosolutions generate". 

The usefulness of the approaches can be underlined by some of the 
results obtained: for example, negative values of the autovalues indicate that 
the finite elements are improper and do not respect any fundamental 
properties of the structure they should approximate. Again the transverse 
deformability causes complex autovalues. A treatment on the nature of the 
autovalues and on their property results interesting. From a practical point of 
view, once the autosolutions have been found, each structural problem (in a 
linear field) can be resolved with their linear combinations, breaking down 
each system (of forces, for example) applied to a certain section in terms of 
eigenvectors. 

In the case of the approaches to finite elements, it is convenient to refer to a 
module analysed with a stiffness method, searching for the displacement 
states that have the property of the constant 4 relation between the 
displacements and -7 between the external forces applied to the terminals 
(that is, between the stress data at the terminals). 

Let us consider any structure 7 that can be ideally divided into a sequence 
of parts A;, A, ... A, here referred as ‘modula’ such that 

A;0A = P; AVAO ..VAL=A 

and where between 4;.; and A; there is a geometrical similarity relation of 
ratio u, independent from i and a material identity. Each ‘mudule’ 4; 
borders only on its preceding and its following ones. The section of border 
are respectively referred to as ‘front’ and ‘rear’ sections. 

External loads can be applied to a module only on these sections. 

Each part of the whole structure can be analyzed through an approach 
that discretizing the same part introduces a finite number of Lagrangian 
coordinates to describe the displacement of front and rear section. A reduced 
external force can be considered for each of these coordinates. 

Each “module” must meet an identical condition of its displacement and 
force reference systems on the two border sections in order to allow the 
adjoining of “modula” and the definition of displacement and force 
configurations having characteristic propagations. 

The task is now to investigate the condition under which in a structure 
like that considered, it may be possible to find displacement and internal 
force configurations having characteristic propagations, i.e. all displacement 


4 The A. has already discussed the problem with the aim of indicating the feasibility of the Cicala's 
work, in Reference [7].[8]. here largely utilized with amendements and clarification up to the end of the 
chapter. 
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and internal forces at each configuration denoted a reduction by a 
multiplication facto through a “module” of the structure. 


3.3.2 Main definitions 


The configuration of a module is defined by means of Lagrangian 
coordinates (in the displacement space) or multiplication factors ; and of 
the correspondent functions D). 

The functions D; constitute a reference system in the space of the element 
boundary displacements, that are consistent with the nature of the 
introduced analysis approach. In an analogous way the internal forces are 
defined by means of Lagrangian coordinates fi and function F). The 
functions constitute a reference system in the space of the forces. In a 
reference system (comprehensive of displacement and forces) usually the 
scalar product between displacement and force give 


(D, F) ai loci 
fa Cus 


i.e. the scalar product among correspondent functions gives a unity of 
mechanical work. 
Let us define 
u=|u;|= displacement column matrix. 
f = | fi|= external reduced forces applied to the structure 
(column matrix) = /,....,m (degrees of freedom 
of the structure). 
The following relation is valid: 


ku=f. 


In the problem we are interested in, the structure displacements that don’t 
interest the boundary cross section can be eliminated by means of the 
correspondent equilibrium equation, where the external forces are zero 
(homogeneous equation). If this operation is done, the system of 
equilibrium, involving only displacements and forces interesting the two 
boundary cross sections, can be put into the form 


[k;]lw4]=|f| , =... 21), (i=1,..,2m) (3.27) 
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where the stiffness matrix £ can be considered as a column matrix composed 
by row matrices A, Each element of A; is a reduced force per unit 
displacement: 


A; 
k=|4, sapa; (3.27) 
Am 
Obviously the k matrix is symmetric and semipositive definite. 
Rigid body displacements 
Let UT ,U) succes Ul, gi U, be a reference system in the space 


of the module rigid body displacement. There is no loss in generality if we 
suppose that the orthogonality 


uz'uz=0 , (a#B) 


is verified. A rigid body displacement can be expanded as follow: 
ACI (3.28) 
a=l 
The conditions of singularity: 
Aju' =0 (3.29) 


holds independently from the choice of u' . 
The condition implies also that 


detk,|=0 


as it necessary in order to have non zero rigid body displacement without 
external forces. 
The following condition of symmetry 


Ik,.|=[|k,.| (3.30) 


holds as result of the reciprocity theorem. 


3 And also \min k| = 0, where min k is each minor of 4, obtained by nullifying displacements in 


such a manner that the structure is still capable of performing rigid body displacements. 





| 
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The virtual work equation states that external force system f is equilibrated 
if its virtual work on a system of rigid body displacement u” is equal to 
zero, independently from the choice of wu”. 


Therefore the condition for the equilibrium is 
,=| ku pE u=0. (3.31) 


The external forces given by (3.27) constitute an equilibrated system. 
This statement can be demonstrated writing L,, of eq. (3.31) into the form 


Li=u"k'u"=0 
and recalling (3.29) and (3.20). 


Module behaviour 


If system (3.27) is considered with the condition (3.31), due to the rigid 
body degrees of freedom, a static problem is 


ur R, 
Mo! DI : 
pero esere Det DM = MS . (3.32) 
i U,.1 Wjr 
Bledinimari [4] i 


U>m a 


The submatrix [1] of the matrix &, is definite positive (and not singular, of 
course). 


If u' is a rigid body displacement, which satisfies the” assigned 
values, we can put 


u=u"+u°, 
where the unknowns w° represent the deformation when conditions on the 
u" are taken into account. It is of course u” = 0 (j=1,....,1). 
We have 
ku°.= f. 


MODERNITY OF PLACIDO CICALA’S THOUGHT 53 


and the problem can be divided into two independent parts. 
The first part is the determination of the unknown displacements u° by 


means of the system 


[4] coll Soda (3.33) 
"7 da 


ei |:|=|:|. (3.34) 


The set of known w" should make unfeasible whatever rigid body 
displacement. If it isn’t the case it must be possible to have from (3.32) 
displacements different from zero with f=0 (homogeneous system). 

The det[[4]| must be in such a case equal to zero. On the contrary, if 


det[4] #0, the system (3.32) gives w‘’s identically zero only if the f ‘s are 


not identically null. 


3.3.3 Characteristic system 
Eingensolution determination 


We can consider a module as a structure realising a relation between 
displacements and external (to the module) forces belonging only to the 
bounds. Let us now choose on the rear section displacement and force 
reference system in such a manner that they coincide with displacement and 
force system of the front section of the similar subsequent modules. 

If such condition is satisfied, we can indicate with w and f respectively 
displacements and forces on the front section and with w, / the 
corresponding ones on the rear section, and write (1) in the form: 
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n A 


A - LB UE pl (3.35) 
B'iC|ul |f 


A 


where: A is a square symmetric matrix of rank m,C is a square 


gr A 


symmetric matrix of rank m, B is a square matrix of rank m, v and w are 


A 


column matrix (of m elements) fand are column matrix (of m 
elements). 

In order a system derived from (3.35) may result invariant through 
physically similar modula it is necessary to introduce the new variables 


iù 


2 
c 


pr ta dal (pt 
c Cc 


EE 9 e; 
[» 


where c is a length. Through obvious operations system (3.35) becomes: 














A  4B 
CB ci” -|/ (3.36) 
dini 
Ko d 
where 
A A G A 
n A da . Ha 3 LEI 
6 E Cc G 
The matrices 
NR | A_i 4B 
k=|4-+8 and k=|B"T C 
rm! dI Lilia 
B_:C Hi 4 


are singular (3.29), but if we consider a submatrix of 2m-r rank or less, 
where r is the number of rigid body degree of freedom (r < m), it result 
definite positive and, obviously, not singular. Therefore the two submatrix A 
and C are definite positive. A characteristic system has 


u=du , f=-3f 


and the equilibrium equation can be written 
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A B 
pi È u - tà 

2 —|[Aul --4 i 
Pe ff 








This equation system can be reduced, through operation, into the following 
one: 
1 B' 
A+ > tABHK == 
J7, À pl 


u=0. (3.37) 








This system is linear homogenous, as unknowns u are concerned, and has 
non zero solution if and only if it is. 
C 1 B' 
A+—+A4BU1+—— 
4 À 
This is an algebraic equation of order 2m in the eigenvalue À.. 
If u=1, the equation can be written 


det =0. (3.38) 











È +C+4B4+B" {u}=0, (3.39) 
that has solution w # 0, only if 
detl4+C+4B+B" =0. (3.40) 








In the most general case, such equation admits complex solutions. 
If B' = B and 4=/, (8) and (9) reduce to 


{w}=0, (341) 





a+c+(2+1)5 
À 





det =0. (3.42) 





a+ca[2+1)a 
A 





Equation (3.38) (and also (3.40)) is an algebraic equation of 2m-order in the 
unknown 4. In the C-space there are 2m (eigevalues) 4 and to each of them 


corresponds an eigenvector u. Obviously if a complex root À exists, the 


conjugate root A there exists too. 
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In the particular case that B' = B eq. (3.40) has real roots for the variable 
I 
unknown 4 + of Such roots are positive if -B is positive definite (A and C 


are definite positive). Another consideration is to be made, related to the 








properties of the considered structures. In fact, the matrix , related to 
r CENTRE si B 
the case B' = B, results semidefinite positive. The problem C {u}= 0 








admits also solutions which are different from zero and are related to the 
rigid body motion. It follows that, if 4% is the eigenvector relate to the 
eigenvalue 4, 

Ta Ea = 1 7 


“=0 or u (A+C)u+2u Bu>0. 








1 


b 
But (A+) (444) Bu=0, therefore 24122 (if B is 


definite negative). 
Therefore it is: 42 1; the values 4=1 corresponds to rigid body motions 
whose related forces are zero. 


U; fr 
Sec. “Rear” state a) Ur,=zU, 
U, a = 
state b) U,j}=U} 
U,, = 0 
Sec. “Front” 
Uf 1 
Figura 3.5 


To note that not all the the finite elements implementations preserve the 
matrix 8 to be negative definite. In such a case one can have negative values 
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I 
for deg and therefore of 4. A better understanding can be reached by 


considering a box which admits, among his characteriste values, 4 to which 
a displacement vector u,; is associated. If one applies a state a) of 
displacements (see Fig.3.5). 

Let ff. be the correspondent vectors of external forces.. 


Let us apply also a displacement state b) which makes complete the previous 
one related to 4 (see Fig. 3.5): 


The correspondent externals forces are: f,,= 4f,, ,f= — 4, - 


SM 


3 U. 
It must be ( being 4 = — = — 


U, VE VW, ja i A, 


Exception made for the case 4=1, which falls in f,, = f,, leading to an 


indetermined form, when Z#1, the vector f,,e f,, must be different each 


i a - 
other only for a constant factor: f.,=@f,,. with x da 


2, 
lea 





a-A=4A-0},a= 


Therefore if a >0->>0;ifa<0-+ 4<0and vivecersa. 
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In those cases in which 4>0, also 2 >0 (as is the case when 8 is 
symmetric and definite negative). 


To notice that the adopted positive signs for f.,. f,,€ @, the positive 


values of Z mean that in correspondence to “rear”, for displacements which 
are zero, the applied forces have opposite sign with respect to the applied to 
“front” section. Therefore, when ,=0 with f,#0, the wu, 


displacement have the same sign of u,. 


A structure which can be divided in more boxes can fulfill or not such a 
property. It can happen that the approximated model of the structures does 
not fulfill the properties which is fulfilled by the real structures. 

Let consider the box depicted in Fig.3.6, analyzed in the framework of the 
reinforced shell. 


Rear 





(1) (2) (3) 


Figura 3.6 


MODERNITY OF PLACIDO CICALA’S THOUGHT 59 


For sake of semplicity we consider the case P=1, A=1,/=1, 431. 

The unit virtual system a) shows an external virtual work 2u,-2u,, which 
coincides to twice displacement of (2) with respect to (1) at the front 
section. 


"= 


Figura 3.7a) 





Figura 3.7 b) 


By  applying the virtual work equations, one  obtains 
l | ; 

MH; ATO , so that w,—w,>0. On the other hand the unit 
Ù) 
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virtual system b) gives the displacement of 2° for respect tol. One can 
7 ] ] 
obtain. u-u=-— 
i Gs. 6E 


and u, — u; results to be a function of s, (see Fig. 3.8). 





Figura 3.8 


The same results can be obtained by fixing s and varying A. It can be seen 
that the sign of @ depends from the geometry which form the base of the 
matematical model. A more sofisticated mathamatical model as those related 
to “shear lag”, will lead to always positive a. 


Eigensolution properties in the complte field 


To a couple of conjugate complex roots 4;, A, , it corresponds a couple of 


conjugate complex eigenvectors %, , uz u, . The following relation 
holds 
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Ù g 
Fptelir ica 1: ripoiiia u,=0. (3.43) 
12. ring iù 


This proves that u 2, is the eigenvector associated to À, ; 


If one exchanges between them rear and front sections, and says A the 
unknown eigenvalue, eq (35) becomes 
y 


urb: (3.44) 








PIASTRA È Bu+% = 
4 A 


2 


A PO ALL | 
The determinantal equation in A gives a root 4, wr) for each 4, The 


system (36) gives for such 4 the eigenvector u'=0 d, where Cl is a 
multiplication factor. 

When B#B' .wecanput B=B,+B,, B,=Bj , B=-Bî. 
The system 3.37) can now be written: 


Pet Ap + la Bf FALLE. B\{u}=0. (345) 
Au Àu° 





u 








If we consider for each solution Z an associated value , Z that can be 
complex that: 


si i di ò 
Au + = -n+ Mis, Au — Lalla pali =-s", 
À pl Te dl 











i 


the following identity holds: 


det = det 7 = 





A PACAA 2 la + [da Db 
Hu Ap Àp° 





I 
Spioici B|=det H!. 
pi 


62 ETTORE ANTONA 


l 
Hence if 4 is solution of det H=0, — is solution of det 7 '=0. 
À 


When 4 =1, obviously e 
À 


] i 

"a It is necessary to note here that when 
Li 

B# B' .the eigenvector associated to Z is different from the associated 


to A, also in the case 4 =1. This result is easy to be verified, starting from 
$B7k 


When 28° # B e u=1, one obtains 
] l 
A+C+|4+—-|B+[{4-—]B, 
À À 


where B=B,+B,, Bi =B, . Bi =-B,.Inthis case the unknown 


À can be complex. Being 


(If 


the previous equation can be written as 


{u}=0, (3.46) 








2 
a+c+[44+1)p, + (4-3) -2B,|{u}=0, 
À À 
whose determinantal equation can be expressed in terms of the variables 
(4 + 1) ; turning out to the variables 4 one notices that the inverse of each 
root is still a solution of the determinantal equation. Not proper solutions 


could appear by solving the previous equation. In fact, 


(444) —-2 =|- (444) —-2 
À À 
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In reality this fact does not happen and the number of the roots for 


1 
(4 Pi 2) are m while those for 4 are 2m°. 


The roots for 4 , if complex, appear in subsets of 4 roots. This is due to 
the presence of the inverse and of the complex conjugate of each roots. 


: ipa 
However, if the number of the roots (4 + 1) is even, at least one roots 


l 
and at least two values (4 e Fi, must be real. 
In the equation (3.46) where A, C, Bo are symmetric and B; is 


i ] : 
antysymmetric, let us consider z = ai ; we obtain 


a+cs(+1)a (2-4) 
z z 


By introducing the symbol x, which states for both 4 and z, the two matrices 
of (3.46) and (3.47) are sum of a common symmetric part plus an 
antysimmetric part, whose second contribution appear as a transposed 


{u}=0. (3.47) 








I 
matrix. That is È — 1) appear in (3.46), while 
D 


l gi pi 
(» La = È a appear in (3.47). 
* x 
The determinantal equation of the two matrices is therefore the same. 


si ha ] | 
That is for each root 4; it exists the root — , and its complex conjugate. 


td 


On the analysis of the reality conditions 

Let us distinguish the displacements of a section into displacements 
orthogonal to the section itself u, and displacements in the plane of the 
section u,. Let us refer to the case (7) where in the submatrix 8 the elements 


concerning only orthogonal displacements v, and forces (f,) are symmetric. 


See note 1) 
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Referring to the part ii) of the eigenvalue problem, if all the deformation 
displacements, that are involved in elements of 8 that are non-symmetric, 
are zeroed and if the corresponding forces are disregarded because not 
interested in the characteristic propagation, the problem can be reduced to a 
real positive eigenvalue problem. 

In a structure when the hypothesis (i) is verified and when there are 
perfectly rigid frames, the condition of no deformation displacements is 
verified by the u and w’, displacements, that under the made hypothesis are 
displacements to which can be associated elements of B that are 
antisymmetrix, and the module has real /s. 


On the symmetry of sub-matrix B 


If one considers the “frames” to be infinitely rigid in their plane, the 
submatrix B results to be symmetrical; this instead does not occur if one 
considers their deformability, always on their plane. 

The fact that the matrix B is not symmetrical means that, taking into 
consideration its position inside the matrix of stiffness of the module, the 


node forces F induced by the displacements U are different from the node 


forces F induced by the displacements U . 

In other words, given a set of displacements applied to an extremity of the 
module, the node forces induced at the other extremity depend not only on 
the displacements themselves, but also on which of the two extremities they 
have been applied to. 

If one considers a module made up of a single panel (Fig.3.9), and 
consider the problem of obtaining the characteristic displacements at its 
extremes, taking any two opposite sides as terminal sections, that special 
submatrix that we have identified as B does not result to be symmetrical. 

The dissymmetry of the matrix B is not caused by the arrangement of the 
panels in the module, but was already implicit in the stiffness matrix of a 
single element. 

If one fixes the directions ( and not the sides) of the displacements for 
each node, all the elements of the stiffness matrix are determined at least by 
the sign. Having fixed the direction, (understood as the straight line of 
action) to a node displacement, taking a side as positive or negative in fact 
means placing a different sign from the elements of the stiffness matrix that 
are correlated to the displacement or the corresponding force. Let's in fact 
change the side of the k-th displacement; this means to change all the 
elements of the k-th line of the stiffness matrix ( the positive side of the 
force in fact changes) and all the elements of the k-th column ( in that the 





MODERNITY OF PLACIDO CICALA’S THOUGHT 65 


side of the displacement changes, and if a positive displacement first gave 
certain effects, a displacement that is positive now gives opposite effects to 
the side). Incidentally, the element of the main diagonal remains unchanged, 
undergoing two changes in sign. 

Let’s consider the module made up of a single panel; to say that the 
submatrix B is not symmetrical means, in more concrete terms, that the 
displacements of the points of an extreme section have, as node forces, 
repercussions on the other section depending on which of the extreme 
section is considered. Let's now suppose to have adopted the same reference 
system at the two extremities of the module and that the panel that makes up 
the module is symmetric to a transversal plane, that is, by reversing the 
module and exchanging the two extremities, a panel that is identical to the 
previous one is obtained ( this is logical if the panel is made up of a single 
rectangular finite element, while it is no longer true in the case where it is 
made up of the joining of two triangular finite elements). 


E Îi 





b DE 
I 
Lins 
Fi 
Fig. 3.9 a), b), c) 

Let°s take the panel represented in the Fig.3.9a as an example. Let's 
consider only two degrees of freedom for each node; the positive sides of 
the forces and of the displacements are those indicated. 

We now want to see what would be the effects of two corresponding 
positive displacements, at the two extremities of the panel. Let°s consider 
the displacements | and I and let's suppose that the node forces respectvely 
inl and | are indicated in the Figs.3.9b,3.9c with the dashed vectors. 
Fig.3.9c is obtainable from Fig.3.9b by turning it by 180° and inverting all 
the sides of the designed vectors. 
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We can now see how, at a component level, referring to the adopted 
reference system, the two cases give different results. 
As far as the first case is concerned, it is possible to state that the 





displacements generate two forces F and F, , both of which are negative to 
the reference system; (-a) and (-b) being their components. 

In the second case, having adopted a still positive displacement 
according to convention, has led to, with respects to the modules, a situation 


that one can call “reversed”, that is, if the displacement U, tends to do 


compression, the displacement 10) now tends to do tension. As already 


mentioned, we can in fact pass from one situation to the other by reversing 
the module and by inverting the sides of the vectors. 


At component level with the positive displacement U; the two 
components of the node forces are no longer both negative, but, always 
according to the references adopted, /) is positive while F} is negative. If 


in both cases the same value or the displacement has been adopted, for 
example unitary, is it right to expect that the modules of the corresponding 


forces are, in the two cases, equal (that is, that the module of the F is 


equal to that of the F, and also for the other two forces). In the case in 


which the module presents symmetry to a parallel axis at its extremity side, 
which is the equivalent of saying that, by reversing the module one obtains 
an identical module to the preceding one. 

In the case in which this occurs, for example when the panel is made up 
of a rectangular finite element, it is therefore right to expect that the 
modules of the corresponding forces are equal, and therefore, in the second 


case, one obtains (a) and (b) as components of the node forces F and F, 


according to the adopted reference. As a conclusion, it is possible to say that 
the components of the corresponding forces are different in the two cases: 
and this is equivalent to saying that the matrix B is not symmetrical ( in one 
case in fact the node forces were found through the matrix B, and in the 
other through its transposed matrix, and we have found different results). 
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3.4 Properties of the eigensolutions and analysis of the characteristes 
systems 


The characteristic systems have interesting properties that allow their 
utilization in various aspect of the structural analysis. For sake of brevity we 
shall limit our esposition to the summarize of some of such properties. On 
eigensolution is concerned with rigid body displacements it is demonstrated 
that among the 2 m configuration having characteristic propagation ther 
exists a 2 r dimension space of solutions, deriving from rigid body 
displacement, without external forces, [8]. In fact the rigid body 
displacements the determinantal equations gives the solution 


NEO 
À 


with an r— multiplicity in A, and X = 1 with a 2r — multiplicity, such roots 
give u'=u, and f = f'=0. Eigenvalues that are not associated to rigid 
body displacements or other force-free displacements are real, if B = B° 


[8]. 


If B # B' the eigenvalues can be complex. To real values 





À 

- #] there 
are associated force distributions whose force and moment resultants on 
each section are null. In fact if such resultants were non null, they owing the 


4 uz albronng Je 
respect — # 1, should create a non equilibrium situation [8]. 


There is a complete separation between free rigid body and deformable 
body eigensolutions. In particular deformable body eigensolutions don’t 


work on each section displacements of the free rigid body type, [8]. 
(a) (b) 


ti 


Let us consider the eigenvectors and w", correspondent 


) 1 
respectively to À, e —. 
A set of m-r eigenvectors, different one from another, and picked-up 


from the whole set 


1) (a) b (b) 
Vi rice avan 


9°m-r? a°m-r? 
is composed of linearly independent vectors [8]. Obviously an analogous 
property there exists as the eigenvectors f‘’ and Fl are concerned. It is 
(a) 


2/ m-r 


easy to conclude that the eigenvectors f‘,..... can be used to expand 


as linear combinations each self-equilibrated force distribution introduced 
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(5) (5) 


on the first section. The eigenvectors f;'”,....., f;,, can be used to expand 


as linear combination each self-equilibrated force distribution introduced on 
the other section, [8]. 


Two eigenvectors u, and u, , if 4, # lt in the space of scalar product 
Ax 
ui By u;, are orthogonal when B;=0: ui Bou;=0. 
The demostration follows from the Betti theorem, [8]. 


Two eigenvectors %, and %,, when to = Lesa are orthogonal in the space 


H da 


si - 
of scalar product u, Bru: u,B,u=0. 


4. CALCULUS OF VARIATIONS 


4.1 Introduction 


Cicala’s approach to the calculus of variations stands out for an 
originality of planning and a great incidence of personal contributions, 
particularly on the indirect methods, and puts at the engineer’s disposal the 
instruments to find maxima or minima of integral expressions and the 
unknown functions that generate such extrema, without the need of notions 
of superior analysis, with particular reference to the cases of only one 
independent variable. 

As it is known, in its foundations the problem received various 
contributions. At a certain point the approaches divided themselves into 
“direct methods” and “indirect methods” and there emerged the need of 
including among the solutions those constituted of the junctions of 
etherogeneus arcs and of considering problems with variable terminal points 
too. The necessity to present the complex corpus of knowledge and methods 
necessary to solve the problems of calculus of variations and optimal control 
in a coherent vision and without the need of superior analysis is, we can 
sustain it, at the basis of Pontrjiagin’s proposal to adopt his “principle” too. 
Cicala’s formulation puts one in a position of total awareness about the 
meaning and the implications of the mathematical conditions to be satisfied. 
Various authors spoke about Cicala’s contribution: for example Kopp, [12], 
observes that the maximum principle derives directly from an extension of 
Weierstrass” and Cicala’s works. 
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A brief review, coherent with our purposes on origins and developments 

of the variations calculus must do reference to various points of view: 

e the various and very frequent enlargements in the complexity of the 
considered problems, 

e the refinement process of test criteria or necessary and sufficient 
conditions for the existence of maxima (minima) among solutions of 
differential problems. 


4.1.1 Evolution of the considered problems 


Newton proposed the problem of the solid of revolution of least drag 
(1686). The quantity to be minimized is the integral depending on an 
unknown function: 


"de 7 
fa CE de (4.1) 
a (TH? 


with a and d given constants. Eq. (4.1) is coherent with a Newtonian theory 
of the fluidynamic drag. The solution obtained by Newton himself was not 
so correct and exaustive, even if starting points can be accepted. 
Weierstrass” analysis (late nineteenth century) provided a reason why such 
solution was not satisfactory from a theoretical point of view. 

The problem proposed and solved by John Bernoulli, said 
“brachistochrone problem” (1696): “find the curve to be followed by a 
heavy body to descend, without friction, from a point A to a point 8 in least 
time”, asks to find a function y(x) which yields the minimum value of the 


integral 


{ Lod a (42) 
V 


The solution obtained by Bernoulli wasn't rigorous and induced discussions 
among the contemporary scientists. In particular James Bernoulli was able 
to find a way of solution applicable to a more general class of problems. 
These problems caused the research of general solutions, not bound by the 





An aero-fluidynamies problem. 


70 ETTORE ANTONA 


specific problem; form these the so-called classical theory of the calculus of 
variations was born. 

James Bernoulli (1697) proposed the “isoperimetric problem”, in which the 
function y(x), apart from being a minimizer for a functional, should satisfy 


the additional integral condition 


fg.) = cost, (4.3) 


so enlarging the field of the variational calculus. 

The isoperimetric problem can be considered a particular case of the so 
called “Lagrange’s problem”, where the unknown function must satisfy a 
differential equation. If we define terminals a and b, the problem considered 
in so far has fixed terminals. The evolution of the nature of the problems 
considered also not fixed terminals. Another enlargement was the case of 
pieced solutions. 

As the complexity of the problem is regarded, the evolution considered 
successively two or more functions y,, two or more conditions @,, with 


terminals subjected to conditions or fixed*. Problems considering functions f 
and @,, including higher y, derivatives, only apparently are of enlarged 


complexity. In fact, introducing subsidiary functions y,;, they can be 
reduced to formulations considering only first y; derivatives. Further, 


Lagrange’s problem can be considered a particular case of Mayers problem. 
Mayer’s formulation of the maxima (minima) problem can be considered the 
more general and permits to pose any problem under the parametric form as 
the unknown functions are considered. 

Hamilton’s formulation is of particular interest, upon which the methods 
of optimal control are based, a discipline developed relatively recently, 
which finds a great motivation in the problems of today’s science and 
technique”. 


4.1.2 Evolution of the test criteria 
Considering as an example the problem of determining the function y(t) 


for which an integral 


8 î i 3 3 ; 
Other kind of problems here not considered are those with two or more independent variables. 

° No mention is made here of the Dirichlet®s problem, starting point of the functional calculus, because it 

considers two independent variables and such kind of problems is not considered by Cicala. 
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E ff (1, y,V)di (4.4) 


is maximum (or minimum), the so called “Euler’s test” was the first to be 
introduced: 


ni bt pe Rf) (4.5) 


It derives obviously from the condition dg = 0 and was obtained in its 


final meaning of a necessary condition through the solution of the 
brachistocrone problem due to John Bernoulli, based on geometrical 
considerations not exempt from intuitive components, and the extension due 
to Euler. 

Later Lagrange, by means of the variation concept, gave a rigorous 
approach to the problem, taking it again in the differential calculus field. 
The distinetion between extremals and maxima (or minima) was made clear, 
in particular by Ruesch. Legendre indicated the necessity of considering also 


$° g in order to demonstrate the existence of maxima (or minima) and gave 


a brilliant demonstration of a condition on 0° f/4° , which must have a 


definite sign. An important observation of Lagrange made clear that also 
such Legendre’s condition is necessary but not sufficient, because it isn't 
always possible to operate a (mathematical) transformation that conduces 


from é°g tod f/3° . This particular mathematical problem was found by 


Jacobi, who deduced a new necessary condition, deriving from Sg, for a 


certain time regarded erroneously as sufficient if in conjunetion with the 
Legendre test. Jacobi’ s condition requires that: 

1. the extremal supposed to give a maximum (or minimum), 

2. the envelope of the extremal having the same origin of 1, 
do not have contact points in the interval a-b. 

In a third period of time, due to the contribution of Weierstrass, Darboux 
and Kneser, other important clarifications were obtained. In particular, in 
1879 Weierstrass demonstrated that Jacobi’s condition was not sufficient, 
by means of a very deep criticism on the adequacy of the way of thinking 
used in the analysis of maxima and minima of the differential calculus. In 
particular he observed that the presence of } in the expression of f, makes 
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not sufficient the sign definiteness of 6° g. Even if small perturbations of y 
are considered, strong variations of y can be involved. 


He introduced a Weierstrass function £, and, as Tonelli said ([11], Vol.1, 
pag. 18), reached the conclusion that for an extremal 


1. Legendre’s condition, Li >0 (<0), 


2. Jacobi’s condition, 
3. Weierstrass” condition E£=> 0 (<0) 


in the whole constitute a sufficient condition for being a maximum (or 
minimum).This is not applicable from the most general point of view, 
because only a class of strong variations were considered. 

In fact: in ref. [10] one can find the statement that the before mentioned 
limitations don’t permit the necessary condition of Weierstrass to be a 
sufficient one. Our opinion is that Weierstrass” condition on the variations 
taken into account regards not only the problem of Legendre - Jacobi, 
concerning the second variation, but also the problem of the first variation 
which drives to Euler-Lagrange’s equations. Weierstrass adopted, as a basis 
for its analysis, functions satisfying Euler - Lagrange’s equations, or in other 
words extremals having 6g = 0, in problems with fixed terminals, if weak 
variations are applied. Considering the class of strong variations that we can 
call “jogged lines”, he proposed the test E>0  (<). If satisfied, such test 
ensures: 

1. that dg = 0 also for jogged variations, 


2. that the extremal is a maximum (minimum). 

Therefore under Jacobi”s condition Weierstrass” test include Legendre’s 
test, but not Jacobi”s one. 

Tonelli.[11], with a very important contribution, dated at the beginning 
of this century, gave a demonstration of the general validity of Weierstrass” 
criterion, taking advantage of the property of semicontinuity of the 
funetionals to be maximized (or minimized), in correspondence of the 
maximizing (or minimizing) function f. 

Also the introduction of the concept of “field of extremals” and the 
consideration of a parametric form for the involved functions and integrals 
are due to Weierstrass. The extremals of a “field” do have the same starting 
point and different terminal point. Both “fields of extremals” and 
“parametric form” are largely referred to in Cicala’s approach. 

An exaustive presentation of the argument should distinguish between i) 
absolute and ii) relative extrema, even if the condition of ii) must be 
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satisfied, as necessary conditions, also by i). Here the conditions for relative 
maxima (or minima) are taken into consideration. 

Other concepts not considered here are those related to “weak” extrema, 
distinguished from the “strong ones”. They are not necessaries, in 
consideration of our aims. 

Other authors and other points of view are not listed here, even if 
interesting and important, because not interested in Cicala”s approach. 

To show Legendre’s necessary condition, let the quantity to be 
minimized be that defined by: 


g= (f@y,3)d (i=1...,n); 4.6) 


then, if we call Y;(Y) the coordinates of the solution of the problem, a 
necessary condition along all the solution is: 


df(1,Y,Y) 


sn > 0 Vi e[t',1°|]. (4.7) 


It is evident from 4.7) that it springs from a condition of the type È g>0, 


in strict analogy with the differential calculus. 

To show Jacobi®s necessary condition, let us consider the space of the 
solution y; , and examine the case of a variational problem for which the 
starting point A is fixed but not the endpoint Z. Obvoiusly the solution has a 
degree of freedom and the whole of the solutions (different for a constant C) 
is called extremal field. If moreover all the solutions are non intersecting, at 


least in the interval £# — #‘ of interest, we are in front of a central extremal 
field or extremal star. Jacobi”s condition requires that the solution candidate 
to be minimizing is part of an extremal star. Stated in more useful terms for 
the applications, the enunciation is different (but in all equivalent to that 
given here): let us consider an extremal field with its origin in A; the 
trajectories of this field could have an envelope or not, that is a curve 
pertaining to all the trajectories and tangent to them. If there exists, Jacobi”s 
condition says that the minimizing trajectory must not have common points 
with the envelope, exception made of eventually the point A (at least in the 
interval 17 — 1‘). We therefore can say that if the extremal under exam 
touches the envelope, then certainly it does not minimize the quantity g. 

The demonstration of Legendre’s and Jacobi’s tests require large 
amounts of mathematical deductions, here not necessary. It is sufficient in 
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consideration of our purposes to summarize some concepts, that in order to 
facilitate the exposition, we will refer to the case of only one function y as in 
(4.4), without integral conditions (see [10], chapter 1). 

The quantity to be maximized is given by 4.4), where y(f) is the 


- è a A B PESO 
extremal which must respects given terminal values y° and y”.A variation 


y=y+%v, where v' = v" = 0 must give the first variation 6g = 0. If we 


want to take into account the behaviour of the second variations, we can 
operate a variation of the first variation y : 


y =V+e7=yY+&+€N, 


where y+ 7) is a varied y, satisfying as y the condition 4.5), and not 


satisfying the values y° and y at the terminals. The first variation of a 
first variation is a second variation, and y = y + €77 must be the solution of 


a problem having the same y° but a different y” (in other words we can 
say “an extremal of the star having A as apex”). By means of mathematical 
operations it is possible to demonstrate that the second variation is 
proportional to 





If the parenthesis is positive (its square surely is positive), second 


9 


variation is defined positive if always A >0 in all A-B (test of 





Legendre). The condition (i- và) > 0 is verified if the parenthesis is 
) 


not zero. The value 77° is obviously null; if 77 can be chosen proportional 
to v, the parenthesis vanishes and Legendre”s test loses its meaning. But this 
is impossible if Jacobi’s test is verified: “considering the envelope of the 
extremal star with apex in A, the extremal y must have no tangent points 
within A-B, terminal B included”. 

We can leave to the reader the extension to the case of several y,. 
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4.2. Cicala’s approach 


As a basis for his contributions Cicala adopted Mayer’s formulation, 
which considers the problem, representative of a great number of cases, of 
finding a curve such that: 

e alongita given ordinary differential system is satisfied; 

e the starting and ending points A , Z satisfy a group of assigned 
conditions; 

e the quantity g which depends upon the collocation of A , Z, if calculated 

through it, reaches a minimum value. 

We are considering a differential problem in which all the quantities 
depend upon only one variable, say /, indicated also with “j=0”. The 
derivative is indicated by mean of }, is. 


The problem can be formulated in the following way: let Y; be the 


coordinates of the solution curve in a (n+/)-space which comprise the points 
A, Z. The differential system is written as 


Pa(;:);)=0; j=0.l..n, i=l.n a=l.m, (4.8) 


m equations of type 4.8), which contain general functions of the n 
coordinates y,, of the n derivatives y, and of time (y, ort). The 


conditions for the points at the extrema of the solution come from 


v;(91.37)=0 Broli.ps (4.9) 
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and the quantity to be minimized is 
8=g9/:37): (4.10) 


FIG. 4.1 





FIG. 4.2 


The variable yo must be monotonically increasing or decreasing, otherwise 
the solution sholud be found joining arcs where / satisfies this condition: 
moreover we don’t consider limits on the quantities y, : they are free to 
assume any value is necessary to follow the optimal trajectory. In the system 
4.8) we can have non derived variables; these, which often present 
themselves as control variables, will be called x, and are in the number of /. 
Let us consider first the problem deriving from the first variations; then the 
system 4.8) will take the form: 


Pa(V;:V;X,)=0; j=0,l..n, i=l.n, a=l.m+l, k=1.l 
(4.11) 


Lagrange’s problem and the isoperimetric problem are included in this 
formulation. The problem is a differential one because we have no more 
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integrals and it is possible to solve it writing equations which permit one to 
determine all the unknowns present in (4.8), (4.9) and (4.10). 

The method used is that of Lagrange’s multipliers, applied opportunely. 
It is possible to transform the problem, where we have differential 
equations, into a bounded ordinary minimum problem. In an original 
approach particularly apt for an engineering use, Cicala proceeds with 


discretizing the equations (4.8), dividing the interval tt" ina number N 
of equal intervals. Doing so the differential system is substituted by a system 
of mN algebraic equations which, together with (4.9), form a set of 
equations of bound for the function g. We then do a limit operation, letting 
N go to infinity. The main results can be summarized as follows: 


Z; =dF/d; 
z:=dF/d,, 
i (4.12) 
Z0 7 ZiVi» 
F=A,@g; 


IAAD E (4.13) 


Id = Di 
0G|1Ò4î =-25, (4.14) 
G=g+A;Yg. 


The 4.12) are Eulero-Lagrange’s equations, in another form; it is 
sufficient to combine the first with the second to obtain them. The 4.13) was 
obtained with further manipulations and represents a result of very great 
conceptual value. It is called “incremental relation” and lets one obtain a 
link between the total variations, denoted with the small delta symbol, of the 
values at the estrema passing from an extremal to a neighboring curve. 
Furthermore it has a role in calculus of variations similar to the equation of 
virtual works in analytical mechanics. The equations 4.14) are called scaling 
relations, because they let one fix the initial value of the multipliers and, 
then, close the solution of the problem by choosing a particular extremal 
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between all of them obtained from the solution of the system 4.8), 4.9) and 
4.12) and defined without the initial values of the multipliers. 
If moreover one among the 4.9) is in the form: 


re N) SOL (4.9?) 


stating that the point A must be on a given ipersurface, then if we apply 
4.14) we obtain the so-called /ransversality relation which states that the 
solution starts perpendicularly to the surface. If for example we have: 


P = f0,y)=0, (4.97) 
we then obtain: 
z(F/9,)=z(d/d). (4.15) 


In the case of parametric representation of the locus of end point, for 
instance a curve in a three dimensional space: 


P,=vf()=0, 
P.=Y.-fiW)=0, (4.16) 
P,=y:-fi()=0, 


the transversality condition has the form: 


z;:f+2,f,+2;f;=0. (4.17) 


This holds particularly if we have pieced solutions, topic introduced by 
Cicala in a general shape, to which he gave fundamental contributions. 

The number a of the independent differential equations must be a<n: if 
a=n we are in presence of a normal differential system and no optimization 
is possible (and only n conditions at the terminal could be done). 

On the other hand the number p of the conditions 4.9), must be p < 2n. 


If p=2n+1 the optimum problem has only one degree of freedom and the 


choice of the optimum is to be performed among arcs that verify equations 
4.8) but are not extremals, exception made for the solution itself. If 
p<2n+1 it is possible to construct the star of the extremals having the 


same terminal A and all the terminals Z that are permitted by equations 
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(4.9). The system of stars starting from all the permitted A, generates a field 
of extremals. AIl the extremals of a field verify conditions (4.8) and (4.9) 
and have first variations equal to zero, as 2n+/ terminal values are 
considered fixed. 

Weierstrass”!° necessary condition is introduced by Cicala in an original 
manner. We now consider a problem formulated according to Mayer, of the 
type defined by eq. 4.8), where the quantity to be estremized can be posed in 
the form: 


girini (4.18) 


where / is a function of the relevant coordinates; these are practically those 
variables for which the z, in 4.12); are not identically zero. 
The scaling condition is: 


A 
AE st (4.19) 


toni PAN 
Ii 





where y, is such that Cile does not vanish in the considered 
h h 
region. 
With reference to fig. (4.1), we consider an extremal path (i.e. a path 
verifying Euler-Lagrange’s tests, that, it is useful to remember, were 
obtained considering only weak variations of the y,) which starts from a 


certain point A and ends in a point Z in a problem with fixed terminals. The 
path C is not extremal and arrives at the same terminal Z: it is relative to a 
variation of the extremal and starts from a point that is different from A only 


» . . A z : 
as the value / is considered (fixed terminals concern y; and y; , exception 


made for the values of % and we have just imposed the value of h° ). The 
path C can include weak and also, more important, strong variations: if the 
extremal path is a maximum (minimum), the projection of C on the plane /-t 
doesn’t intercept the projection of the extremal itself except at Z. For a point 
Y of C we can consider the extremum with fixed terminal that starts from A 
(value of 4 not fixed) and arrive at Y. If the Jacobi test is satisfied, such 
extremal doesn't intersect the A-Z one. 


Lo; dra soul CE0E 
If joined to other conditions it can result as sufficient. 
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Let us denote with an asterisk the quantities referring to this last small 
are (arbitrary, but verifying the field equation 4.8). 


So we denote by dv, = dt, dvi = },di° the components of the arc Y- 


PA 
Then, according to Weierstrass the following relation holds: 


ég = Edi”, (4.20) 
where £ is Weierstrass’ excess function: 
E=z;(),-3}). (4.21) 


The 4.20) is interpreted as follows: the variation of the quantity g passing 
from the “right” optimal path A-Z to the alternative path B-Y-Z depends on 


the sign of E and of dt° and on the time during which the path follows the 
arc of generic shape carrying from Y to Z. 

Weierstrass” condition springs from these considerations: if we want the 
extremal under examination is a minimizer for g, it is necessary that along 
such extremal it is always £>0,if it is verified the complementary condition 


dt° > 0, that is any path the solution can follow, differing from the arc A-Z, 
the quantity g must continuosly receive positive contributions greater than 
those it could have had, at the same instant, if it had followed the optimal 
arc. 

This condition, introduced by Cicala in such an original way and of great 
conceptual and pratical relevance, springs from autonomous reasoning 
without the need of the theory of ordinary maxima and minima; so we can 
see in Weierstrass an early attempt to face some problems of calculus of 
variations starting, with a different point of view, from advanced results as 
Euler’s equations and the extremal field. 


The condition di° > 0 gives a new and relevant importance to Jacobi”s 
test. In fact, in order to registrate a change of sign of dt° in A-Z, it is 
necessary a crossing of B-Y with A-Z. The Jacobi test prevent from the 
intersection in A-Z of two extremals. 

It is useful to note that the definite positivity of Edt°  permits the 
conclusion that g is semicontinuous (in Tonelli?s meaning) at the extremal 
considered. So both the requirements deriving from ée and 6°g are 
satisfied, also for strong variations, even if the extremals considered are 
obtained with weak variations. 
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The fields of pieced extremals are very important for the practical 
applications. It is proved that even in the points of discontinuity, the paths in 
a field of pieced extremals most always have tangents forming angles 


a <90° with respect to the preceeding ones (the curve goes on). 

On the other hand Weierstrass” condition £>0 reduces to Legendre’s one 
in the case of weak variations. We propose here a demonstration of such a 
statement in a simple case of two functions y and g of an independent 
variable x, without non-derivated variables, and using for Weierstrass” 
condition the way of thinking of Cicala. 

With reference to figure (4.1), if only weak variations are taken into 


account, the line having differences x° — x and y° — y infinitesimal of the 


first order and can be considered another extremal D-Z. The previously 
considered extremal A-Z can be considered a variation of D-Z. The points of 
A-Z are defined by the values x, y, y, g, £ . The points of D-Z are defined by 


the values x°,y°,y°,g",é .A point X on A-Z can be the starting point of 
an extremal E-X, where E has the values of x + of D. The value of hf is 


obviously free (for respect to h? ) in order to permit to E-X of being an 
extremal. 


Applying the condition of Cicala, we can write for kh: 


* * . . 
where z, and z, are derived from a function 


* 


F=# -flx.5), 


valid on D-Z, variation of the function 
Fi = E- f(x,Y,)) , 


valid on A-Z. 
The condition F° = F =0 give: 
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Lx -n)-4(5-))=0. (4.22) 


#-9-La'- 0-1 


n 


It is obviously 


The value of z; =dF" P. d° can be obtained by means of an expansion 


starting from F: 





sdtevoig@i gli X°G lennotgargiiogi b 
z; (uit) e peptldirzd): 


Weierstrass’ excess function can be written: 








E=ale-#)-Z6-3Y Ta at 
i (4.23) 


The extremal E-X can be chosen, without any loose in generality, in such a 
manner that x° = x. Egs. 4.22) and 4.23), in such a condition give: 





DI 


We have £ > 0 if it is = > 0, which is Legendre’s condition for the 





simple proposed problem. The extensions to the general case having i 
greater than two and the presence of non derived variables x, (note that 


Euler?s conditions give F/R, = 0) is left to the reader. 


4.3 Problems with pieced solutions 


The search for the extremals for problems with pieced solutions, studied 
systematically by Cicala, is of great importance for the engineering 
applications; the solving curve of the problem is constituted of ares of 
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estremals each satisfying a given differential set and boundary conditions, 
which have common extrema. The sources of discontinuity can be of 
different nature. For example the solution often isn’t free to reach all the 
possible values, but must satisfy one or more conditions of the type @, 2 0, 


which we will call g/oba! inequalities. On the other hand, there could be 
discontinuities in the formulation of the problem (equations with different 
definitions if one or more variable exceed a given value). Moreover we can 
have constraints on some quantities (it is of particular interest the case 
where the constraints are on the non-derivated quantities, which often 
represent controls). 

If we denote with a,b,...,z the indexes which identify the various ares of 
extremals, we have generally speaking z different sets of equations of the 
kind of 4.11) and then of Euler-Lagrange equations of kind (4.12). At the 
junction between two subsequent arcs the equations (4.13) will assume the 
form: 


y Za Ab Za 
dG1 di +01 97 = 20 200, (4.24) 
where the first apex denotes the point under consideration (initial, A or final, 
Z) on the arc denoted by the second apex (a=first arc, b=second are and so 
on...). 

As usual all the quantities to be minimized appear in G together with a// the 
possible constraints present at the nodes. The equation (4.24) contains cases 
remarkable and well known; if the coordinates y, at the passage from one arc 
to another are left free (we don’t fix the transition point) then we obtain 
Erdmann-Weierstrass” condition 


Az; =0, (4.25) 


while if the coordinates vit are subject to a constraint ‘P, > 0 we have 


another well known result: 
OP. 
Az; = Ul, () } (4.26) 
Ab 


The 4.24), 4.25) and 4.26) are known as corner conditions. The former 
(Euler-Lagrange’s equations) are necessary and pose conditions to the 
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solution in the form of differential equations to be satisfied. There exist 
other necessary conditions, among which we have Legendre’s one and 
Jacobi”s one. They are based on the analysis of the second variation of the 
functional to be minimized, because it was a persuasion of their authors (but 
more generally it was a common thinking at that epoch among the 
mathematicians) that the problems of calculus of variations could be solved 
in analogy with differential calculus. Therefore, because necessary 
conditions exist for the maxima and minima of functions (based on the 
analysis of the first derivative) and sufficient conditions (based on the 
analysis of the sign of the second derivative) it was common to think that 
this way would have been similarly fruitful of results. The reality was very 
different and the first man who understood this was Weierstrass who 
formulated by proceeding in a totally different way, the necessary condition 
which carries his name . Practically the second order necessary conditions 
which we will examine allow us, among all the solutions solving Mayer's 
problem, as defined by 4.9), 4.10), 4.11) , 4.12), to eliminate those not 
satisfying these additional conditions. 


4.4. Constrained quantities as conditions 


As to Lagrange’s problem with discontinuities, Cicala shows how, 
through Mayer’s formulation, a treatment is possible even in the presence of 
constrained quantities (they can originate a discontinuity in the solution). 
The starting point is the same: determine a solution y,(f) minimizing the 


integral 


[ea . (4.6rip) 
If we consider 
) sn pb. Va ini tv a t), (4.27) 


and write, according to Mayer, 
F=A[f01:009-y.]+4,01-A. (4.28) 


g=yi vi, (4.29) 
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we obtain: 


-_ 
Ill 
= 
5 


= Va t0911), (4.30)a pe 


and if we eliminate ff} and A, we have the well known Euler-Lagrange’s 
equation: 


VECI 


î . (4.31) 
di, di 


Now we focus our attention on the f function: for assigned values of £ and y, 
it is a function of f}. The function (8) will play a primary role in all that 
follows and will be recalled in the following with the term of f-curve. Each 
point of the plane /,y; (that is the trajectory plane) identifies one f-curve in 
the plane f, f. The direction of the tangent to an extremal y, = f(1) in a 
point P of it (i.e. the quantity y, |, ) identifies one point on the f-curve which 
will be called indicative point. So, according to the equations 31).pc the 
value 4, represents the slope of the f-curve in the indicative point. Next we 
will see that A, plays a fundamental role in the process of drawing the field 


of pieced extremal paths. We note that, if B is the value of y, for a path in 


the plane / , y; of the trajectories alternative respect to the optimal one but 
always passing through P, it characterize on the f-curve a different 


indicative point, whose coordinates are (B', i "). Then, the excess function 
locally is given by: 


E=A,05-V)+A,0-V)=f -S-4P - BI). (4.32) 


If we put the condition £>0, we obtain the important geometrie condition 
that if we want the solution of the problem to be a minimizer, the whole f- 
curve must lie above the tangent to itself in the indicative point. 

We study a very common case if we establish some limit to the possible 
variation of /, that is the velocity y,. We will write 8, < 2 < 58, where /, 


and 8, can be function of / and y,. The solution will be composed of ares 
where the condition 2 = 8, holds, ares where 8 = f, and ares where we 


haven’t any bound on /. Where the first limit is active we will write 
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F=A(f-.)+A,4-B.)+A;(B-B;), (4.33) 


where the bound was introduced through the multiplier A3. Being 8 costant, 
the solution y, can be calculated with a direct integration of the equation 


Y;= f;(t,y,). We then obtain y,from »,=f; Euler-Lagrange’s 
equations are given by: 





ott Pi Ri quei 
Ar=a 4a A35+4,, 4,30, (4.34) be 
and eliminating À we obtain 


À, + Pusendi, (4.35) 
db, di 


from which we can calculate the multipliers, particularly 4. The equation 
4.34), shows that a/ong the constrained arc (with the indicative point fixed 
at one estremum of the f-curve), the local tangent is different from the value 
of À;. At a Junction between a constrained arc and a free one it is possible to 
show, by using the corner conditions, various properties that for sake of 
brevity aren’t recalled here. 
AII what we said about the are 8 = 8, holds true for the are 8 = £,, if we 
change opportunely the symbols. 

It is not useless at this point to resume how can we do the entire process 
of drawing the extremal path: given the function f, 8} e #8 then, draw the f- 
curve and find the tangents to the extrema of the permitted arc for the values 
of , here called s,(f,y,) and s,(f,y,), and calculate the values of 4; 


according to 4.35); then follow the rule of table I. 
TABLE I 


iii, Pal é bb 
if-4><f; + p=B, 
if, >sS, > B=P, 
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An original contribution of great relevance due to Cicala is the 
interpretation, which we can may define “graphic-engineeristic” of 
Weierstrass” condition. The excess function is posed in the form: 


E=z($,-5))20. (4.36) 


in which the )) are quantities subject to verify, together with the x;, the 


conditions: 
Pa (VI) = 0, (4.37) 


while all the quantities calculated with the values of x and vi are marked 


with an asterisk. 

By a genial coupling of the possibilities offered from this representation and 
from the concept of extremal field, Cicala shows the possibility of a 
graphical discussion of the most hidden problems of optimization. 


4.4. Second variations and their analysis 


In Mayer’s formulation the quantity to be minimized is a function g (see 
eq. 4.18) of values in A and Z, conditioned by equations ‘P, = 0 (see eq. 


4.9) of values in A and Z. This allows us to formulate an optimization 
problem as a maximum (minimum) problem, starting from a function G (see 
eq. 4.14) containing g an the P, ’s. Let M be an extremal that maximize (or 


minimize) g with the conditions P, =0, and N an extremal in the 
neighbourhood of M, verifying the conditions ‘P, = 0. If we consider 


quantities 77, subjected to the conditions 


(34) =0 (#8=1,...,p) 
3, n; = ulo42 


in order M to be a maximum (minimum) the condition: 
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9 


(É°)am <o >0) (,k=1,...,n) 
Fi 





must be respected, or in other words the definite negativity (positivity) of a 
quadratic form must be demonstrated. 

Another way of considering the same problem is to introduce parameters 
£, (r=1,..,5) for each of the degrees of freedom, concerning the 


location of the end points (the simplest problem has only one degree of 
freedom, that of the quantity to be maximized or minimized). If 5° is the 


optimal configuration, 


Og c 
be be EE, = Er EE, 
r k 





must be of definite sign. At this point Cicala introduces an impressive 
suggestion that permits avoiding the serious difficulties that one could meet 
in the before mentioned way. Such suggestion does a brilliant use of the fact 
that the first variation at a varied point is the second variation in the non 
varied point, exception made for the remainder. The analysis is so reduced 
to the &/de, in varied points. The problem is solved by the following 


steps, where s is the number of degrees of freedom: 


I. g,£,€, must have definite sign; this is verified if: 
2. the determinant | La | has definite sign, with all the determinants obtained 


cancelling in each of them the last line and the last row; this is verified if: 
3, 1/e, (x/0e,) the definite sign requested in 1) and 2), for the solutions 


E 41 “a E 4) sila Es E 0 


satisfying the conditions | 3 
&/0e,=@/den =...= &/0, =0 


As far as the analysis of the second variations is concerned, Cicala 
introduces a brilliant way of thinking, where the optimum path is considered 
as belonging to a family of extremals satisfying egs. 4.8) and 4.9). Each path 
of the family is defined by a value of a parameter b and has coordinates 
expressed in the form y;(f,b), and 6 being independent along the arc. At 
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the terminal, obviously, { as well as the y, are functions of d, and we can 


write 
Dal r Dis) 
LL =(y/ D e+ T=A,Z), 
(È (0) +) pp): 1342 


where y=d4/®=n;, );=d;/a, di/db= rt. 
With F=4,@, and F= {Fa, 


The demonstration of such result is here omitted. It makes use of the 
properties @, = F = GF/0%,=0, of the integration by parts, of Euler's 


; dF 
equations and disregard the summands as Fr", F7*, bognd: dt and 
so on, because they are null as effects of 4.8) and of Euler’s conditions. It is 
dF 
also db =0 and consequentially is zero the sum of the second member, 


but are zero in it the various summands. We can write: 


vi (E) La 


Eq. 4.8) and Euler’s equations for the intermediate ares hold for each 
extremal. On the other hand the parameter d generates a family of extremal 
arcs. Therefore it is @,=0 F=0 d7/01,=0 d/d;-î;=0 for 


all the values of b. Operating we obtain: 





I 
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dF 2 dy; 7 dy, i 

—=401|, SARE iL TORO 

db db bl, 
SCAN Ar i 
db db\' JT ab db|,) 








At the extremals A and Z is @,, = P, = 0 forall values of d. Consequently 








we have: 
d dI, dy; 
Pa 14 Pa ip 
dy, dj, dy; db 
Therefore we have: 
de_ GB bj GG Dj 
dh dî db d°' db’ 

















db? db db db 


3 db 3 db 
J di 


Concluding, we obtain: 


di 





+ ini 
db 


#) 


db db 


d°G_ d FALLE d È 223 


E ii Ape 
Buon 





dy° 
db 


F dG 
Considering that a =0, because of 4.8), and d =0, because of the 


condition on the first order variation of G, and considering also that the 
second member holds for any choice of the derivatives dy È [db, we can 


conclude that: 





(4.38) 
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Considering the second derivative, one obtains: 








di) db dlbi |) db’ 


ch ngi; Lr ent DE 
ti di 
(4.39) 


In this demonstration, here reported in a way easier to be understood, one 
can appreciate Cicala’s sharpness of mind, in particular for not having 
2 4 To) PERTIII! î 
4 in the dF /db expression, in order to dispose separately the 





zeroed z,); 
summands of z;(dy, /db) and to demonstrate 4.38) and 4.39). 


As the problem involves r degree of freedom, d "n db? can be reduced to 
a homogeneous quadratic form in ” parameters. The positiveness of such 
quadratic form furnishes the minimum condition for the extremal under 
consideration. We can now add some arguments to Cicala”s demonstration, 
that furnishes sufficient conditions in order to verify that among the totality 


of the extremals coherent with (9), the considered extremal y; generates a 


minimum for g. 
IDA 0 i STARS 
The more general variation of y; can be decomposed in a variation of 


another extremal, having terminals coherent with 4.9), and a variation 
coherent with 4.8) having fixed terminals. This last variation generates 
second variations definitely positive if Weierstrass” test and Jacobi’s test are 
fulfilled. In conclusion: 


1. Definite positivity of d° g/db° ; 
2. Weierstrass test, 
3. Jacobi test, 


considered globally, constitute a sufficient condition for Di to be a 


minimum. 

The formulation of the analysis of second variations in the extremal field, 
due to Cicala, allow us to obtain also in another way the necessary Jacobi”s 
conditions for fixed end points problems. In fact the limit situation to avoid 
the presence of a coniugate point of the terminal A in the arc A-Z is revealed 
by solutions 77; of the necessary problem vanishing at both end points. 
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4.6. A sample problem 


Cicala’s treatment of the problem in general is indeed imposing and rich 
in sharp observations. Since we haven’t the aim of discussing the problem in 
a general way, we will refer to a particularly simple example. 

The following case will also permit the comparison between Cicala’s 
method and that which uses the “maximum principle”. 

Let us consider the problem of designing a beam of minimum weigth whose 
free torsional swinging must have a minimum assigned frequence. If we can 
apply St. Venant’s theory, the equazions defining the problem are: 


M=R9, 
M+Jo?9=0, (4.40) 
W=w, 


being A the torsional stiffness, M the torque, ./ the moment of inertia about 
the elastic axis of the mass of a unitary length, @ the pulsation, w the weigth 
of the unit length. Let the derivatives be calculated respect to the axial 
coordinate £. W is the weigth of the length until the section 4. The section 
t=0 is built-in while the other end is free. The quantity to be minimized is 
then: 


g=W°, 
and the function F is: 


F=A,(w-W)+4,(M/R-9)+4,(M+J0?9). 


From this we can derive Euler-Lagrange”s equation: 


À, =0, 
À,+J0°4,=0, (4.41) 
A,=4,/R. 


' 
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The dimension of the elements of the sections are supposed known as 
functions of / and of one variable x (e.g. the thickness of the wall in a thin 
walled section). So w and R are known functions of x and /. Let it be the 
effect of x on J negligible, while the last is a known function of £. Moreover 
it is prescribed the condition R > R,, where the right hand term is an 


m? 
assigned function of . The solution, generally speaking, will be composed 
of free ares where R> R,, and constrained ares where R= R,. For the 


m? 


former type the variation of F respect to x gives: 
w'=A,MR/ R°, (4.42) 


which is useful to determine the values of x along the beam, while for the 
latter x is given by the equation of the constrain. 

We can use Weierstrass° condition to build the extremal, which here is 
written as follows: 





3a 
Krk (4.43) 


and it can be visualized with the help of the graph of the curve 
w= f(1/R). 

The condition requires that this latter is all above the straight line whose 
inclination is -9 and passing through the indicative point. While it is 
R> R,,.-q is the local inclination of the diagram. On the contrary along the 
constrained arcs -qg is greater than the local inclination of the diagram at the 
extremum R= R,. 

At the junction of the two arcs the following physical conditions hold: 
A9= AM = AM=0 and then the corner condition give: 


Aw+qgA(1/R)=0. (4.44) 


We have to specifiy the boundary conditions for the differential problem 
given by 4.40) and 4.41). 
As to the former system, we prescribe the physical conditions 


M° = 9'=0. The transversality conditions for the multipliers give 
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A,=1e 4f=A43=0. We only mention the fact that we have the same 
formal differential system holding for both the state variables and the 4°s. 
This fact carries the result that along the beam is 9 / A, = constant, thus 
facilitating enormously the solution of the problem. 


It is interesting to show how the same problem can be faced using the 
“maximum principle”. 
The status equation are the same field equations 4.40) (the variables that 
appear with their derivatives are called the status variables) with the same 
quantity g to be minimized; from these we form the Hamiltonian in the 
following manner: 


H=Y A,%,=A,w+A,M!/R-A,J0°9 
j=i 


and, from this, we formulate the adjoint equations in the form: 


À=-V.H, 


which here is written: 


(we note the different sign respect to the 4.41); the transversality conditions 
are the following: 


A =(dol&)f , (4.45) 


where the derivatives are to be calculated respect to those state variables 
whose terminal value isn’t known “a priori”. If instead one variable is 
specified, generally speaking, the corresponding multiplier, for {= #7, will 
be non zero. The 4.45) here give: 
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2 È 
A. 
AJ 


Moreover, if for f = 1" we have not prescribed the value of some state 
variable, here the corresponding multiplier is zero. In the case under exam 


this gives 45 = 0. 
If we have free arcs it is necessary that the derivative of the Hamiltonian 
respect to the non-derivated variable x (said “the controller”) is zero: 


si =w'-A,MR/ R°=0 


and so we have obtained the 4.42). 

To discuss about the constrained ares we have to call the “maximum 
principle”, which states it is necessary for the optimality of a certain 
solution that it is a minimizer of the Hamiltonian respect to all the other 
possible solutions compatible with the status equations only. 

We can note here, at least intuitively, the link existing between this 
enunciation and Weierstrass” condition as used in Cicala’s point of view: 
really, Weierstrass” function could be considered conceptually a sort of 
derivative of the Hamiltonian respect to the trajectories (which are governed 
lastly by the controller x): if we move along the optimal trajectory , it is 
allowed to expect that any other will give higher values of g. 

This intuitive speech can be formulated in the following way: let us 
consider the variation of 7 respect to the controller and say it is positive; we 
will obtain the condition: 


= -w+a Mi) > 0, 
RR 


which is equal, apart from the equal sign, with the (4.43); so we find the 
same inequality pilots the composition of the extremal in both cases. 

It is possible, starting always from the “principle”, to show that along a 
constrained are R= R, the tangent to the curve w= f(1/R) in the 


indicative point must differ from the quantity -g. Indeed, the allowed 
variations of the controller x must be so that the corresponding variations in 
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R are necessarily positive, because this latter cannot descend below Rm. So, 
if we assume a monotonic law x > AR, if we want the Hamiltonian to be 
minimized it is necessary that: 


20. 


lesa 
In our case this writes 
w'-A,MRI R°>0 


or 


Lastly, the corner conditions at the junction (for /=1,) are obtained from the 
following equations: 


ASSI Da (Ad x 
At) 





o (4.46) 
H:=Hnoaagf, 
à 


l 


where NM(x(f,),t,)=0is the equation (sometimes a system) which binds 
some state variable in the interior points. 

The preceeding system, when applied to the case under examination gives, 
for the multipliers, the result that the transition occurs without jumps in the 
adjoint variables 4, while the Hamiltonian is continuous there. Because we 
often suppose that the curve w = f(1/R) does not permit jumps in the 
controller (possible, if it has straight arcs, or single isolated points) the 
second equation of (4.46) gives the (4.44). 
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4.7. Modernity of Cicala’s thought. 


Briefly, the modernity of Cicala’s thought about this subject springs from 
his effort to arrive through an approach with great notes of originality at a 
formulation which allows to face the real problems of the engineering, even 
the most articulate ones. 

The applications testify of the potential development, some of great 
relevance, because applied to problems such those of the aerospace activity. 

As we can see, the problem exposed before contains a condition 
R= Rin: a type belived for a long time impossible to be treated with 
Calculus of Variations and for which we should have to use Pontrjiagin’s 
“principle”, considered for this reason more general. 

Conditions of this type, which are a reason, together with others, of the 
“pieced” nature of the solutions, can be generally written: 


0,20 (4.47) 


and here the function F is, according to Cicala and in a more general 
approach than the preceeding, 


F=A,@,+ AQ, (4.48) 


in which 4, is zero or not depending on the fact that the arc is in the field 
@,>0 or on the bound @, =0. With this planning Cicala gives one 


ultimate contribution more to the possibility to obtain the same applicability 
of “Pontrjiagin’s Principle” from Mayer’s planning, through the use of 
Weierstrass” condition and the same Cicala’s geometric view. 

As an example, to cut the length of this treatment, let us consider the 
problem of minimum length between two points in a tridimensional space 
with coordinates y,, and in presence of an obstacle 


Y(;)=0. (4.49) 


The F, if f= y, is the path along the trajectory, is given by 


F=Ay+A;(yd,-1)/2. (4.50) 
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From the definitions of 2; we obtain 

z:=-A); 23-74, (4.51) 
which holds in all the field. The G can be written 

G=t° + uy(y), (4.52) 


where x is zero in the space free from the condition y = 0. Apart from a 
scaling condition A, =1, here given as obvious, Euler”s equations on the 
obstacle give 


D y 
A (4.53) 
pEr) 


that is the arcs on the obstacle are “geodesics”. In the free field the arcs are 
straight. For the junction point between a free arc ‘a’ and an arc ‘b’ on the 
obstacle, we have the condition 


sesta. 





(4.54) 


If we multiply both the sides of 52) with the first addendum of the right 
hand side, we obtain 


dae 1 1 7 (4.55) 
i.e. the ares have the same tangent. 


Trajectories of this kind can touch angular points. Even in this case we can 
use “Cicala’s engineeristic approach”. 
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5 Asymptotic approaches 


5.1. General remarks !! 


The asymptotic approaches, to which P. Cicala dedicated many 
publications (from the sixties), concerning principally the thin walled 
structures, constitute an original and revolutionary issue. It can be applied to 
the searching for the approximate solutions when some of the quantities 
involved in the definition of the problem can be considered small with 
respect to the remaining ones (for instance the thickness in the shell or the 
transversal dimensions in the beams). This approach could constitute the 
basic philosophy of many engeneering subjects and many chapters of the 
mathematical physic. It is opposed to the axiomatic approach that starts with 
specific assumptions (sometime ingenious) imposing constraints to the same 
nature of the solution. P.Cicala gave an enormous contribution in this field, 
building a theory in which all the axiomatic approaches are not only 
particular cases, but they are also comprehended in their limits and 
implications. 

The results obtained from St.Venant, Timoshenko, Kirchhoff, Mindlin 
and others, appear to be overhung, at a logical level, from a vision founded 
over the subsequent principal lines (shell analysis is here particularly 
considered): 

- the three-dimensional problem is considered belonging to a family of 
problems generated by a vanishing thickness parameter; the behaviour 
of each unknown along the same thickness is expanded in series of 
functions that constitute a complete base: the coefficients of each series 
result functions of a coordinate system defined on the midsurface and 
the problem is reduced into two dimensions; 

- theeffectof the derivations with respect to the space coordinates, in the 
order of smalIness of the unknowns is taken into consideration ; 

- the order of smallness of the solutions and of the infinite terms of the 
equations is evaluated before solving the equations itself ; 

- the systemsof infinite equations are subdivided into partial sub-systems 
each of them with an even number of terms (in finite number) and with a 
specific order of smalIness, neglecting the other. 


1! This argument was previously discussed by the Author in [15], where the work of P. Cicala was 
considered as a basis for the discussion. Present paragraph and part of the next one are largely derived 
from [15]. 
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Cicala created for each of these operations, genial and original 
procedures in which some surprising aspects can be found : the possibility 
of evaluating the order of smallness of the unknowns before solving the 
system, the type at least semi-convergent of the series obtainable in certain 
refining processes of the solutions and so on. 

The asymptotic approach can advantageously be used in each problem 
reducible to an analytical formulation that maintains its significance when 
an intrinsic parameter tends to zero. Comparing it to the axiomatic approach, 
in which the semplifications suggested by the smallness of certain 
parameters are introduced heuristically as a priori assumptions, there are the 
subsequent advantages: 

i) an approximating level is not imposed by certain axiom regarding 
the solution; 

ii) a sequence of steps of growing level of approximation is determined 
where the axiomatic model are defined as steps in the sequence; 

iii) each step introduces a connection between loads and boundary 
conditions. 

An elementary algebraic linear system is considered as discussed by Cicala 

himself: 


D' 4,03) PV) = 0 e, je (5.1) 

j=1 

x; = Xò, l (0.2) 

The a;,N,, are constant and dò, isa Jj—_th power of the parameter d . As it 


tends to zero, the system can be divided in subsystems that can be solved in 
sequence. Each subsystem contains a group of unknown and known values 
calculated in the previous step. The expressions indicated with (2) are 
obtained where dé is a 4; —th power of é. 
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Ina specific case with n=5 , the coefficients matrix is as: 
0 do dg Wa @ a,sò 


dj0 —dGx;0 
0 db, @y 0 a,;ò |. (5.3) 


dn An dA; 
andò andò azò ay asd 


a,ò6 agdò az3dò day 4 


Applicating the Kramer rule, the order of smalIness of the factor dj in (2) 
corresponding to N,=1, is indicated: 


6,=D,/D (5-4) 


where D is the order of smallness of the determinant and D;, is the order of 
smallness of the cofactor of the ij element in the matrix. For each j the lower 
value is taken with i varying. The order of smallness of each addendum in 
the system is d,6; . The fundamental terms in each equations are those of 


the lower power of é. In the case of N; = 1 , the subsequent expressions are 
obtained: 
6=0,=00;05=6?;6,=0;=d. (5.5) 


Tab.5.1 is obtained considering in the l.h.s only the foundamental terms and 
the known terms. 


Table 5.1. Scheme N;=1. 
















3° -d13X3, 
3° +a,X1= -223X3, 
ò +a3,x1ò -a43X3, 
ò +a52X2d +a51X1ò +as3x3ò, 


A 
+a30X30" 


Xo =>X0, X]j =X8°. X4 =Xò, Xs =X50, X3 =X;8°. 





102 ETTORE ANTONA 





For N;=1, the order of the fundamental terms is shown in the left column; for N;=1 
the order of the unknowns is shown in the bottom row. 





The same table is valid in the case: NW,<5%, N,j<é6, Ns<é, N3<d?. 
Neglecting the r.h.s. terms, that are as small as 6? or more, each equation 


contains only one further unknown. 
The case N;=/, generates the Table n. 5.2 and the relations: 


Table 5.2. Scheme No=1. 





XxX} =Xx6% X4 =X4ò, Xs =X50, X3 = 3 Xo = 6 





For N;=1, the order of the fundamental terms is shown in the left column; for N3=1 
the order of the unknowns is shown in the bottom row. 





6,=6°, 6,=6,=6, 6,=6;=6. (5.6) 


In the case N) =/ the table n.5.2 is obtained and consequently the relations 
(6): the solving steps are different. When a reduced system is solved, the 
values of the unknowns are inserted in the second members considering 
them as known quantities, for a subsequent step that takes into consideration 
the most important part of the known terms itself. The system related to the 
case N=1 and N,<6%, N,j<6, N;5<6, N3<&?, is written in the form 
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Dx" (IL, KX}+{N},)=0 
n=0 


with 7=ò ? and the summation extended to n=1. If the solution is written 


in the form: 
{x}=L2"Bh: (5.7) 


the solving procedure follows the steps : 


The expression (7) shows that in the case of systems of infinite equations, 
the solution {x} are summations of convergent or at least semiconvergent 
series: this is a very important property. 

The algebraic example previously considered is a simulation of a wider 
problem and will be reconsidered in the subsequent section. 
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Figure 5.1- Geometrical support of the model. 


MODERNITY OF PLACIDO CICALA’S THOUGHT 105 





Figure 5.2- In plane and transverse deformation components. 
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Fdkdidt , 
= 





Figure 5.3- In plane and transverse stress components. 
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Figure 5.4- Exponents of h in the generic set of nine equations constituting the 
infinite system in the two-dimensional formulation of linear shell theory: the symbols 
=, 2= indicate that successive integers enter in subsequent blocks; the symbol =1 
indicates that successive integers enter in the preceding blocks; the symbols -1", 0" 
indicate that the same exponents enter in alternate subsequent blocks. 





Figure 5.5- Polyhedron of points (r, s, n). A= (4, 0, 0), B (0, 4, 0), C (8, 0, 2), D (0, 
8, 2), E (10, 0, 4), F (0, 10,4). 
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Figure 5.6. Classes of solutions. 
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5.2. The application of the asymptotic approach to the shells 


The application of the concepts of the asymptotic approach to the shells 
is necessarily based on a theory (linear) of the shell itself. Let us indicate 
with 

z=%,+ Gi (5.8) 

the position vector of a generic point P, with x, a point on the medium 

surface, expressed in curvilinear coordinates £,,é,, with n the unit vector 


normal to the surface, and $ the transverse coordinate. 


There are the derivatives : 


F=A105, S=Y10, S = 


and the definitions: 


, (5.9a) 








k,=k,<k, myfi=k,xk, h,=%,, h=%,. (5.9) 


The fig.5.1 shows the geometrical support of the considered model. Fig.5.2 
and 5.3 give respectively the scheme of inplane and transverse deformation 
components and the scheme of inplane and transverse stress components. 
The stress resultants s,, sy, s , have contravariant components defined by: 


Sa = Sa a +s5uwk, + Sac! 
S,=Sak,+ sky (5.10, b,c) 
S. = Sag + Spell, +S,M 


u,=U a? 
u,=t-h,, (5.11) 
UIEZIUCII 
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The relations of elasticity are given by : 
ks, +(kî, = vmî,)s, +2k7k,ySa = ks, = Em,, (5.12) 
13k;s, + (ki, — vmî,)s, +2k}K,Sa — Vk}s, = Em, o (5.13) 


2(1+ v)m2, 5a +2k, (K25, +k75, +2k,5, —8,)= Em, (K,<u+k,-i°), 


Pad =. 
AA 
Tr I 
k, <U 


(5.14) 

s, = V(Kis, +k}5, +2k,,54) = Em,,n-u°, (5.15) 
N<4k?s,, +k3,5,, = Em, (k,-U° +î-#°)/21+v), (5.16) 
k}s,, +k3,5,, = Em, (k, <-&° +î-2)/2(1+v). (5.17) 


The infinitesimal equilibrium equations of the element dé,,dé,,dC , 


subjected to volume loads as f, is the equation (5.18). The group of 
equations indicated with / represents the planar elasticity equations that is 
(5.12)-(5.14). The group of equations indicated with N are the transversal 
elasticity equations as (5.15)-(5.17). The relation: 


5+5+50+f=0 (5.18) 
give three scalar equations indicated with £ . The nine unknowns : 
U=(U,,U,;u,), S = (5395-54): T = (543945980) (È) 


are: the three displacement components w,, u», w, the three planar stresses 5, 
Ss, Sab, the three transversal stresses sy, Spe, Sc. 

The nine unknowns contained in the U,S,7 symbols, are expressed in series 
of the Legendre polynomials in the variable 24/4, where A is the thickness: 


G,= (5 1 AIA E -R 14), (5.19) 
(£)=I@i+D/h [fede . (5.20) 


The polynomials are indicated in (5.19) with j an integer from zero to 
infinity. The (5.20) is the application of the legendre operator to a generic 
function . Operating in a function V, the component Vj can be obtained: 
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V= DREREDO, Ò 


The same happens for the derivatives V V if the thickness is constant, 
according to: 


W°),=[4j+2/4], 


jH1 


baro da (5.21) 
(9); = LR/Aj DIS + +DA/AÎ+ DIS) aa - (5.22) 


The product (é’V) is simplified by the recursive expression (5.22) with the 
reduction of the exponent n to zero. In the (5.12)-(5.17) the m2,5 is a second 
degree polynomial in 4. The division for this quantity origins a series of 
power of 6 for the first members in (5.12)-(5.17) while the second members 
are linear. The Legendre operator is applicable to the transformed equations 
(5.12)-(5.17) . The /j and Nj, that contains the Legendre component of the 
unknowns as indicated in fig.5.4. The numbers in the table are the / power 
in the addenda. The points are zeros, the presence of the derivative operator 
is indicated with a point or a prime or an apex. The extension of the table for 
lower and upper index in the unknowns components is indicated: 


Î=s0+5/+f, 
(5) ja /1= 0); /Aj+6)-@);/(4j+2), (5.23) 


©=0 @A12=6) ©h16= 6) (5.24) 


The fig.5.4 shows also the scheme for the equations £;., of the equilibrium 
equations. They are rewritten in (5.23) with the notation indicated, avoiding 
the summation of infinite terms considered in (8); referring to the (5.21). 
These vectorial equations contain linear function in £. 

The scheme (fig.5.4) is repeated along the principal diagonal except in the 
initial part where the reductions (5.24) form a particular mathematical 
structure. 

The initial and strong characterization of each shell analysis approach is 
the equation of the midsurface. Another fundamental choice that 
characterizes the analysis is the reference system of the midsurface. In 
consequence of the approach including the application of Legendre’s 


112 ETTORE ANTONA 


operator, one obtains the scheme of the differential system where only the 
unknown terms are considered. If one substitutes the differentiations by X 
and Y the parallel system containing only the unknown terms is worked out. 
An operation of substitution could reduce the parallel system to an 
homogeneous equation containing only an unknown. Each term of such 
equation is of the type: ah"X"Y'° (a is a constant). The polyhedron of the 
type of fig.5.5 can be formed. If X=0(47), X=0(417), each term is of 
the order : i=n+ar+/8. Such an equation defines a plane in the 


polyhedron space. The substitution of the X,Y corresponds to taking into 
account the influences of the two derivations on the order of smalIness of a 
function. The same concept can be expressed considering the lengths of 
variation in the two coordinate directions taking into account their order of 
smallness. As far as the response to external not homogeneous loads is 
concerned, the lengths of variation of the unknown funetions depend on 
those of the external loads themselves. As far as the integrals of the 
homogeneous equations are concerned, the dependence of the order of 
smallness from the derivations is included in the general properties of such 
integrals: for instance in the case of separation of the variables, integrals can 
be put in the exponential form e’, and the derivation produces effects 
depending on the order of smalIness of 4 . Linearity allows superposition of 
effects. Therefore the known terms can be considered one by one. If the 
equation of the differential system where a known term is considered, is the 
last in the operation of substitution, a plane tangent to the polyedron and 
having the two already defined variation indexes imposed by the known 
term, evidentiates the more important terms. Each couple of values of the 
indexes individuates a class od solutions. For each class the parallel system 
gives the 4 — table where the order of smalIness of the various coefficients 
are reported. By the consideration of an unitary known term separately in 
each equation (such term must respect the values of the indexes in the 
considered class), the — reciprocal table can be obtained, where are 
reported line by line the orders of smalIness of the various unknown. Such 
orders are obtained by means of the Kramer rule, considering only the 
determinants at the numerator of the various unknown. The reciprocal 
matrix is also the transpose of the matrix of the parallel system. A_ loading 
state can be analyzed in summands of different classes and a superposition 
can be made (taking into account any difference in the order of smallness of 
the various summand if necessary) in order to obtain a particular integral of 
the not homogeneous system. As this particular integral is concerned it is 
sufficient that the polyedron is touched in a point. The lines (equations) that 
in the 4 reciprocal have the same numbers, constitute the leading set of a 
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type of solutions. More types of solution can be present in the same class. If 
we consider now the homogeneous problem of a class, we must suppose that 
at least one of the line is a combination of the remaining and must be 
discarded. The solution can be obtained in function of an unknown 
considered as known. The Kramer rule applied to this case gives the same 
number obtained in the reciprocal matrix when an external load was applied 
to the equation here discarded. For generating an integral of the 
homogeneous system, at least two unknown terms in an equation and at least 
two equations must intervene in the same integral with the same order of 
smallness of the unknowns. Therefore only a leading set can generate an 
integral of the homogeneous system. 


5.3. Membrane states as a simple example 


Cicala’s solutions and the fundamental concepts regarding the 
asymptotic approach are presented in the subsequent, particularly 
concerning the definition of the membrane state. As dé tends to zero, the 
shell tends to its medium surface where curvilinear coordinates are defined. 
The definition of an order of smalIness of the derivation (or the effect of the 
derivation on the order of smalIness of the variable) is a basic concept. The 


relation are introduced: 


Ordo A, 01/06, =Y. 

The generic term in the bidimensional differential equations is as follows 
Ah"X°Y"V, 

where A is independent from é , % is the thickness and the unknown V is 


derived a and d times in the coordinate respectively. The order of smallIness 
of the term is in the expressions: 


Voss u=n-aa-bp, 

rd, YisidA. 
The numbers & and fare the “variation indexes” of the considered class of 
solutions. This terminology was adopted by Goldenveizer , but the same 
concept was introduced by Cicala independently. It extend to the differential 
operator the property already examinated for the algebric factors. Another 
extension of the cited concepts is the locus of the representative points of 
the various addenda of the equation obtained (by substitution ) reducing the 
system to one unknown. In the case of an infinite system that locus is an 


open region in the tridimensional space of the coordinates : a,b,n>0. In the 
case under examination there is a triangular lower face delimited by the 


(3:25) 
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points: (0,0,0), (40,0), (0,4.0) that corresponds to @=8=0. For this class of 
solutions , the derivatives with respect to the coordinates don’t alter the 
order of smallness of the variables. The triangular face contains the 
representative points of the same order. The table in fig.5.6 contains the é 
exponents of the coefficients in the equations in the case @=8=0here 
considered. The last line contains the exponents of the unknowns. In each 
equations (orizzontal line) the terms, where the summation of orders of the 
coefficients and of the unknown has the lower values, are considered 
fundamental and are in a circle. This classification of the unknowns 
corresponds to a loading term of é° in the equilibrium equations £y . They 
constitute the well known subsystem of the the membrane second order 
stresses, if reduced to the fundamental ones. The corresponding strains 
origin the displacements U, with the /, system similar to £y of the elliptic 
kind and with positive curvature shell. The subsystem N) gives the U, 
components regarding the rotations of the normal fiber. The contributions of 
the 7, are neglecting becouse smaller than é°. this confirms the Kirchoff 
assumption about the membrane state. The third U, component is the 
elongation of the normal fiber due to Poisson’s effect. The system /, is 
consecutive to the previous one and generates the bending and torsion 
stresses due to the rotations indicated previously and the secondary U, 
effects. By means of E, the transversal stresses 7) are obtained. The five 
triads of equations are simulated with the five algebric equations used 
before as an example. The five unknowns corresponds to the first five triads 
of unknowns of the infinite system in the table reported in fig.5.6. The 
solution of the algebric example in the case of N)=/ are in accordance with 
the classification reported in the last line of fig.5.6. The other values of 
order higher than 2 can be easily obtained. The same system with known 
terms in /) follows the scheme indicated in the example with N;=/; Sy, 7). 
are small as é with respect to the bending stresses $,, characterizing an 
inextensional state. These solutions represent specific strains of the shell 
due to boundary loads and displacements of the respective classes. In the 
case of general boundary conditions there are combinations of solutions of 
the different classes. 

The considered problems are essentially the linear ones where the 
homogeneous and complete equations are involved. In the case of complete 
equations, the known terms are subject of an asymptotic analysis. The same 
operation is applied to the boundary conditions, characterizing the solutions 
of homogeneous equations. Different subdivitions of the system of equations 
are obtained for each component (of known term or boundary conditions). A 
stress state, solution of a structural problem, appears to be the summation of 
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various states each of them characterized by a specific order of smallIness of 
its length of variation, and for this by a proper capability of penetrating the 
structure. Each of those states is associated to the corresponding component 
of the boundary conditions or of the known terms. 

The procedure is applicable to static and dynamic problems and the 
preliminary aspects of its extension to the nonlinear cases are posed by 
Cicala. 


5.4 The asymptotic approach in nonlinear problems 


Cicala treated the problems of rod and elastic beams in non-rectilinear 
initial configuration and showing finite displacements on the basis of the 
concepts of asymptotic approach. In this situation Cicala’s deductions are 
really enlightening. The work by Prof. Cicala can be summarized in the 
subsequent principal steps: 

a) Definition of the initial equations. The strain-displacement equations 
and the local equilibrium equations are deduced from the tridimensional 
theory [18],[19]. The basic assumption concerns the smalIness of the 
strain sustained by the structural element, in order to remain in the 
elastic range. For this reason linear constitutive stress-strain equations 
are adopted. No conditions are posed on the magnitude of the 
displacements and rotations. The strain-displacement equations are 
expressed in terms of elongation and rotation of a principal line over 
which a local dispacement, smaller that the other, is superimposed. This 
introduces some simplifications and separates the contribution of the 
global deformation from the contribution of the local deformation in the 
general expressions. The equations indicated here are considered as a 
whole and no substitution is made in order to express them in terms of 
displacements. 

b) The order of smallness of the variable involved is defined. All the 
variables are connected to the small parameter e that gives the level of 
the maximum tension strain in the structural element. It is possible [18] 
to classify the structures on the basis of the order of that parameter. 

c) The order of each coefficient of the equations considered is established; 
it is possible to evaluate the order of smallness of the unknowns, 
consequently each term in the considered equations has the specific 
order and can be compared with the others. The most important terms 
are collected and the others neglected in order to obtain an approximate 
expression for each equation. The evaluation of the important terms is 
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extended also to the nonlinear ones. If nonlinear terms are of the same 
order as the principal ones, they have to be considered in the 
approximate equation. 


This approach give the possibility of taking into account all the necessary 
terms to build a consistent approximate model. Arbitrary simplifications are 
avoided because of the comparing of the order of the different terms. The 
system is more tractable and it doesn’t oblige the operator to manage 
cumbersome equations. The approximate model is consistent with the 
conditions introduced and gives also an idea of the field of applicability of 
the approximation. Different approximate solutions are obtained changing 
the initial conditions as shown in the subsequent section. An important 
concept, utilized in this part of Cicala’s work, is the length of variation 
substantially previously introduced by Goldenveizer. Supposing that W” is 
the derivative of the variable W with respect to the coordinate x , and 
indicating with Z the characteristic length of the structural element, it can be 
written: 


W'=WIL,, L,=Le, 


y71 


where L, is the length of variation of the variable W for that class of 
solutions. The class of solution is characterized by the exponent z. The 
value 4=0 is assumed, concerning the case where the structural element is 
considered to be far from the constraints and from the points of loading 
application, and with smooth variation of the section characteristics. 

Considering now a rod contained in the plane of coordinate a,b having a 
definite initial curvature and with #/L=£. This is a clamped free rod and two 
external loads Fa, Fb, are applied at the free end. The loads are orientated 
towards the two coordinate axes and Ca,Cd are the two projections of the 
rod over the axes. The initial curvature has finite entity as the two lengths 
defined before. Using the expression of the longitudinal strain, the entity of 
the curvature variation is defined in order to have a maximum stress tension 
of order £. The bending moment remains of order * and the two external 
loads are consequently of the same order. The resultant longitudinal load 
obtained by integration of the longitudinal stress over the section, remains of 
order e* and consequently e) = ? if the axes are centroidal. This value is 
negligible if compared to the bending strain that remains of order one. If the 
same rod is initially straight , the curvature variation remains of order 0 and 
the deflection introduces a bending component of the same order 4 as the 
deflection is of order 0. In this situation the elongation ey is negligible. The 
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presence of a new load Fc perpendicular to the plane of the other two, 
introduces a torsion moment. The stresses due to torsion have to be of the 
same order 1 as the bending stresses. For this reason the torsion variation is 
of order 0 as the bending curvature variation and the torsion moment 
remains of order 4 as the other two moments. For the beam configuration 
there are analogous considerations obviously changing the order of the 
various quantities involved. The various cases are summarised in tab.5.3: 





+=straight or small initial curvature 

“= due to bending in orthogonal plane 
If in the beam case, there is the presence of a distributed torsion moment of 
adequate order, the order of the torsion variation become the same as the 
curvature variation that is 0.5. In order to present an appreciative example, 
let us considered here an initially straight beam. The strain-displacement 
equations can be derived in complete form by the general theory of the 
elastic beams [18,19] . Their expressions reduced to the linear terms in the 
unknowns are: 


2E0 = 2(€, +Uoo + A0,6, - A®0,$,). 


(5.26) 


2E,,=U,3+U; 


Here ep,4@,(a=0,1,2) are the elongation and the generalised curvature 
variations of the reference line of the structure respectively; £;;, are the 


strain components and w, are the local displacement components . The 
number after a comma is an abbreviation for the derivative with respect to 
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the same coordinate. Substituting the strain-displacement relations (5.26) 
into the constitutive relations, six stress-displacement relations are obtained: 


Co (V/A+V)(0g +0; +0) = 2G(e, +55 +40,6, - A0}É,), 
O -(V/1+M)(0% +0,;1+0,,) = 2GU, 3 


0), -(V/(1+V)(Cy +07; +0,,)= 2Gu,,, 
22 ( 00 Il 22 PI (5.27) 
Coi = G(U, +0 —- A@96,), 


Oy 3 G(Un, +Uz0 +A@6,), 
0), =G(U,,+U,;), 


where v is the Poisson coefficient and G is the shear modulus of the 
material. Instead of the last equation, the following equation is introduced: 


011,22 +01 20312 > V(#,1+* lo, +0, +2G(Uoo | (5.28). 


Where the * is for the square bracket expression. The equilibrium equations 
are: 


T00,0 + To11 +03 = 0. 
To,0 +11 +91,3 = 0, (5.29) 
T02,0 + 91,1 +9 30. 


Considering that the unit volume load is of higher order than the tension 
involved, the nine equations resulting from (5.27),(5.28) and (5.29), in the 
nine unknowns, permit the resolution of the problem. The choice of the type 
of structure to be investigated is a starting point : let us analyse the beam 
structures, where the parameter are : %/L=g%. The structures here 
considered have a massive section or a closed form section, in which the two 
transversal dimensions are of the same order. Therefore if the elastic 
characteristics are of the same order, the bending stiffness is of the same 
order as the torsional stiffness . In order to avoid the fractional exponent, all 
the orders of smallness are counted twice. For this reason the case into 


consideration has an 4/L=£! value of the classification parameter. On the 
basis of the indicated basic assumptions, the order of the variations of 
curvature of the structure is the same as the thickness . This gives an order 
of £° for the resultant bending moments acting over the section ad gives an 
order £ for the bending stress o. The torsion moment is according to a 
distributed external load gy of an order that assigns to the torsion moment 
the same order 5. Regarding the order of the resultant axial load, it can be 
said that if it is of the order 5 as the other quantities, the elongation of the 





MODERNITY OF PLACIDO CICALA’S THOUGHT 119 


reference line assumes the order £*. In this case the tensile stress is of 
higher order with respect of the bending one and can be neglected in the 
derivation of the model. The resultant behaviour is typically inextensional . 
With these assumptions, the coefficients of the solving equations originate 
the table 5.4, where the derivative with respect to the transversal coordinate 
introduces the same order as the coordinate in the order of the coefficient, 
while the derivative with respect to the longitudinal coordinate introduces an 
order 0 in the coefficient. This is representative of the behaviour of the 
structure far from the loading points, the boundary, and with smooth 
variation of the section dimensions. 

Table5.4 





















































The four last columns in Tab.5.4 represent the order of smallness of the 
coefficients of the variables assumed as input data in the computation. They 
are the elongation of the reference line ey , the bending curvature and the 
twist. If 40, = e?“ = e/, ep = 3, the known term vector becomes: 


DI Si IRTRILRSTI . 


Combining it with the inverse of the order of smallness, the resultant 
classification for the unknowns of the problem is calculated and reported in 
the line signed with ** in the table. Some general considerations can be 
done in order to arrive to the approximate solving system: 

a) the longitudinal strain is reduced to the bending terms and the 
longitudinal stress is such that: 


Gg EA, -A@,6;): (5.30) 


b) the torsion problem is reduced to the St.Venant approximation: 
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Oy = GU A®0$3), Co = G(Un3 +A@06,): Cori +T0,3 #0, (5.31) 


with the usual boundary conditions, 

c) the other stresses can be neglected because of an higher order, 

d) the shear effect is neglected because the corresponding stresses are of an 
order higher than the bending stresses, they are reduced to the torsional 
stresses. 


According to the cited references, the resultant loads and moments are 
calculated referring to centroidal transverse axis. 

The equilibrium equations of an infinitesimal slice of the beam are as 
follows 


Nyo N;®@,+N,@,+ py=0, 

No = N@®,+N,@, + p = 0, 

NrgN;0,,+ N, + Py=0; 

M,°+M,@,-M,0,+9;=(1+€)N,, 

M.,0-M,@+M,@+9,=-(1+€)N,, 

Mio -M,0,+M,0,+9 =0, 
Introducing the previous assumptions, the resultant loads and moments 
(5.32) are such that : v, = £° ; M,='; and the equations are reduced to : 


(5.32) 


Nio+PD,=0, N+p,=0, Nine ia, 
| 0.0 * Po 10 T Pi 20T Pa (5.33) 


M,0+9,=N,, M,,+q9,=N,, Mx0+9=0. 


The approximation consists of the Eulero-Bernoulli case. With the same 
classification of the variables, the elongation ey can reaches the order £ 
The equations are modified as follows: 


2E0 = 2(e, +A@,6, - A@,6,), 
2€0, = Ugi 7 A®06,, (5.34) 
20, =Ug3 +A®p6,- 
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The resultant loads and moments become: w,='.N,,=e°,M,=e'. The 
resultant equilibrium equations are reduced to: 

Noo + Po = 0. 

N,o+No®,+p,=0, 

No -No®, + p, = 0, 


N95 
M,0+9,=N,, ( ) 
M,0+9,=-N, 
Mo0+90=0, 


with the bending moments of the same order as the torsional moment. The 
relations show now a coupling between the longitudinal load and the 
curvature in the transverse equilibrium and the presence of a nonlinearity of 
the Von Karman type. Expressing the elongation and the curvatures of the 
reference line with respect to the three displacements w,v,w and a rotation @, 
in the initial triad, the two conditions showed give the classical linear beam 
model and the classical nonlinear beam model. In fact the input variables 
are: 


e =u,+1/2(V)+w7)+..., 


=D +... 
i (5.36) 


DI SVI... 

If the transversal displacements have an order of smallness higher than the 
section dimensions, that is: u= £Ì ; v,w= e; $=£° the result is: e, =}: 
®, =? and the relations are those characteristic of the linear model. It 


is important to note that the inplane displacement is of one order higher than 
the transverse one. If the transverse displacements are of the same order as 


. . . bp) 
the transverse dimension of the beam, that is: u= 7, v,w= e, $=e£°;the 
result is : e, =&°. @,=£.The approximation corresponds to a nonlinear 


model. 
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5.5. Peculiarity of Cicala’s thought 


As “asymptotic approaches” are presented in the literature methods that 
intend to explain the mathematical foundation to the well known axiomatic 
formulations regarding the two-dimensional and mono-dimensional 
problems. These methods, comprehensive of the linear and nonlinear 
formulations (in the strain-displacements relations), doesn’t usually consider 
the loading points, the constraints, the boundary... 

Cicala’s approach affords successfully the complete problem in the linear 
case, also the accidental circumstances, giving an autonomous way in order 
to reduce the three-dimensional problem to a two-dimensional (for shells) or 
mono-dimensional (for beams) one. In the case of the nonlinear mono- 
dimensional problems, the initial shape and strain of the reference line are 
taken into consideration. 

The actuality of Cicala’s thought is evident: 

a) every simplified model is obtained from a logical process that 
guarantees the de-coupling of the different behaviours neglecting the 
terms of higher order; 

b) the different levels of approximation of each model consider all the 
terms of the same importance. 


The axiomatic approaches are found as particular cases. In the linear 
configuration each problem consists of a summation of different addenda 
founded over a specific class of solutions. Each class of solution is 
characterised by a specific ability to propagate the boundary effects into the 
structure. This is a point of potential developing. A clear vision of the result 
obtained from the asymptotic approach is a useful basis for application of 
discrete methods principally in the coherence with the specific class of 
solution. This improves the ability of the analysis with respect to the 
approximations of local effects reducing the computation time elsewhere. 
Cicala’s approach could be of great utility also in fields of recent interest as 
the multi-layered structures and the investigation of the behaviour between 
the layers in a general theory that considers certain solutions not dependent 
from aprioristic assumptions. 
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Non linearity in concrete structures 


Franco LEVI 


1. Introductory considerations 


This report is meant as a selective abstract of a research study on non 
linear design of concrete structures conducted in recent years with P. 
MARRO and G. VIARA, in two distinct stages: presentation of a working 
document in form of a Bulletin (No. 227) of Eurointernational Committee 
on reinforced concrete (C.E.B.) and a subsequent updated edition of this 
document, which appeared as a sequence of three memoirs published by the 
Academy of Science of Turin [1] [2]. For obvious reasons of space, we shall 
focus our attention on the most salient results, as given in the later version. 


2. Non linear design and plasticity theory 


It is deemed useful to define from the outset the respective positions of 
the two approaches. As it is known, the static theorem of plasticity theory 
demonstrates that bearing capacity of a statically admissible configuration is 
lower than or equal to the limit load. According to the kinematic theorem, 
instead, bearing capacity of a kinematically admissible configuration is 
higher than or equal to limit load. Such results can be justified by intuitive 
reasoning: it is logical for the static method to take on the character of a 
“lower bound” procedure, since it interprets the problem merely in terms of 
equilibrium and neglects to take into account the statically indeterminate 
restraints. Similarly, it is not surprising that the kinematic theorem has an 
“upper bound” nature, as it sets no limits to deformability. In this 
perspective, non linear design, which takes simultaneously into account both 
restraints and deformations, appears as an intermediate approach, closer to 
physical reality. 


Our simplistic interpretation, however, leaves out non secondary aspects. 
In static approach - such as, for instance, the classical “Ties and struts” 
theory which reduces the resisting mechanism to a statically determinate 


126 FRANCO LEVI 


system of ties and struts - it is indispensable to apply reduction factors to 
strength of the struts in order to allow for the inevitable presence of parasitic 
influences - flexural effects and multi-axial states of stress. This heightens 
the lower bound nature of this approach. In kinematic approach, instead, it 
would be necessary to set limits on deformations, both in relation of failure 
and in relation of serviceability conditions. The resulting decrease in bearing 
capacity would reduce the upper bound character of the method and, in this 
case, would bring the results closer to those of non-linear design. It is 
apparent, however, that effects of redundant restraints and deformability 
limits are beyond the working field of plastic theory, hence the need to 
define analytical methods able to ensure simultaneous respect of all 
mechanical conditions at play. 


3. Further prerequisites for a correct non-linear analysis 


The foregoing considerations fail to take into account a number of 
decisive aspects. We are thinking of effects of probabilistic and stochastic 
scattering of data and complications arising from collaboration between two 
different materials, typical of reinforced concrete. 


It is therefore logical to illustrate the devices adopted to overcome these 
difficulties within the framework of the modern “limit state theory”. 


The European code on reinforced concrete (Eurocode 2) adopts the so- 
called “semi-probabilistic approach” which consists of applying partial 
safety factors to strengths and actions. Such factors are higher for concrete 
and variable loads. When applied to statically determinate structures, this 
method defines bearing capacity by comparing design load-effects (arising 
from increased actions) with design strengths (reduced by corrective 
factors). Failure takes place in the section identified as critical one by 
equilibrium equations. In this framework, geometric stochastic variability 
(variation in mechanical properties along the axis) is generally neglected. 
Stochastic variability of loads over time, instead, is taken into account by 
means of simplified devices. 


More complex is the case of mono or multi-axial redundant structures, in 
which local mechanical response affects the distribution of load effects. 
Hence the need to adopt models which take into account, to some extent, the 
unavoidable influence of probabilistic and stochastic factors. The solution 
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adopted in Eurocode 2 consists of attributing initially to all structural 
elements their mean properties. At points where, with increasing loads, 
critical situations - yielding of steel, excessive shortening of concrete - are 
reached, the (reduced) design properties are applied. This assumption, by 
which critical sections are made to coincide with reduced mechanical 
properties, provides implicitly a solution to the problem of geometric 
stochastic variability. The different variabilities of steel and concrete 
properties are taken into account by diversifying the corresponding safety 
coefficients. 


However, this procedure is exposed to several objections. 
a) The reduction of mechanical properties in critical zones, as envisaged by 
EC2, might seem a highly conservative assumption. In actual fact, the 
resulting redistributions do not necessarily conform to real ones and may 
even work against safety. The Eurocode model is therefore questionable. 
b) The discontinuities corresponding to the shift from mean to design 
mechanical properties entail substantial operating inaccuracies, especially in 
case of failure on concrete side. 
c) If, in a real structure, the strengths in critical zones exceed the values 
assumed by the design model, the accompanying load-effects (shear, torsion, 
anchorage forces) might be underestimated. A drawback related to 
geometrical stocasticity. 


To remedy these difficulties, there are two solutions based on advanced 
theoretical procedures. One consists of applying the so-called Level 2 
probabilistic approach, which, by working in hyperspace of probabilistic 
variables, makes it possible to provide a quantitative evaluation of safety 
margin. Unfortunately, the method is not easy to apply and its results are 
affected by deterministic parameters of the problem. Its application to 
current structures would therefore require extensive collections of tables. In 
actual practice, level 2 proves useful solely as a calibration tool in well 
defined conditions. The alternative is based on using numerical programs 
which take into account the variability of mechanical properties of 
individual finite elements. An approach which can be adopted only for 
important structures. 


Recently a relatively simple method has been proposed, to counter some 
of the objections raised to EC2. This is the so-called “mean values Yr” 
method. In practical terms, the initial idea is to reduce concrete strength a 
priori, by attributing to both materials weighted mean strengths with which 
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non-linear analysis is performed. A posteriori, a single safety factor (Yr) is 
applied on strength side, which is the same regardless of type of failure. 
This offers several advantages: 

a) It eliminates the pre-determined coincidence: lower strength — critical 
zones; 

b) It does away with discontinuities in response diagrams (e.g. moments- 
rotations); 

c) It also avoids the risk of underestimating aecompanying actions. Non- 
linear analysis, in fact, is conducted down to usual level: Y. qx + Y gx and 
ensuing load-effects, no longer limited by local presence of design strengths, 
are compared with the strengths previously reduced by application of Yr 
factor. 

Yet the mean values Yr approach is still at the calibration stage. 


4. Problems concerning ultimate limit states 
4.1 Premature yielding 


The three-span T beam illustrated in fig. 1, dimensioned on the basis of 
the elastic moment diagram for the loading conditions shown in the figure, 
and supporting to increasing loads, was subjected to variable transfers of 
moments from the intermediate supports to lateral spans. Bearing capacities 
of beams were calculated with elastic method and with non-linear EC2 
procedure. The results, i.e. bearing capacity and safety factor for a service 
load of 15.4 KN/m, are collected in Table 1. 
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Table 1 - Three span T beam. Influence of redistribution at the design stage 


redistrib. from the bearing capacities global safety 
support to the lat. span |(kN/m) related to rare 


FERZ=I= Msint NL 
Ms 
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21.93 
22.07 


In Table 2, instead, we give the values of the stresses acting in concrete and 
steel in critical sections under service load, as well as values of the 
parameters defined below. 





Table 2 - Three span T beam - Influence of redistribution at the design stage 
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In Table 2, the symbols PH and C refer to ultimate limit state and stand for 
plastic hinge and collapse, respectively. è, ò,, and è, parameters represent, 
respectively: 


ò : extentof design redistribution 


. ultimate available resisting moment 


ò 


u 


ultimate elastic moment 
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. serviceability moment with crack induced redistribution 


ò 


$ 


elastic moment in serviceability conditions 


The content of the tables shows that: 


Table 1 

- Linear bearing capacity decreases in accordance with the redistribution 
performed 

- The non linear capacity, instead, remains virtually unaltered. 

(Hence the concomitant variation of safety factors) 


Table 2 

The steel stress in serviceability conditions at support section increases 
with increasing design redistribution (reduction of the è factor and hence of 
parameter ò,). However, already for è = 0.80, the stress is  virtually 
equivalent to design yield limit (fya = 430/1.15 N/mm). By reducing ò even 
further, the yield limit is therefore exceeded and the elastic design 
underlying Table 2 ceases to be valid. The example, however, shows the 
great redistribution capability of load effects at ultimate limit state 
characterising non-linear response. It should be noted however that 
redistributions entail concomitant increases of steel stresses at serviceability 
level. In extreme cases, the design yield limit might be passed even before 
rare load combination is reached. In actual practice, this may happen 
whenever the loading conditions differ markedly from the conditions 
assumed as basis of elastic dimensioning. 


4.2 Aspects of response at failure stage 


Analysis of a large number of examples has identified two different 
forms of behaviour at ultimate limit state. If, at collapse, either one single 
hinge is formed or several hinges that do not interact, and therefore do not 
create a pseudo-mechanism, we find that, in last stage of response, load may 
still increase. This is due not so much to the weak evolution of moment 
acting where collapse is localised, but to outstanding large redistribution. 
We shall say that, in these circumstances, ultimate response is pseudo- 
elastic. In this case, bearing capacity is relatively sensitive to the amplitude 
of available plastic rotation. Similarly, the presence of an imposed 
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deformation (settlement, heat effect) generally has non-negligible effects on 
bearing capacity. 


The second possibility is the presence of closely spaced hinges that tend 
to create a mechanism. In this case, in the last stage of response, the load 
increases very slowly and a reduction in rotation capacity of the hinges will 
have modest effects. Conversely, the strains in plastified zones are 
noticeable and the effects of imposed deformations on bearing capacity 
become negligible. 


Pseudo-elastic response is typical of beams with a small number of spans 
or structures in which load is concentrated in restricted zones. Pseudo- 
plastic response, instead, is observed in the presence of large-sized 
structures subjected to distributed loads. 


4.3 Approximate non-linear checking procedures at ultimate limit state 


Eurocode 2 suggests two approaches: the linear method with 
redistribution (LR) and the plastic method (P). The former permits a 
redistribution of load-effects at ultimate limit state governed by the 
following formula: 


8>044+125Ì: 820.7 
d (1) 


where x/d defines the position of neutral axis at ultimate limit state, after 
redistribution. EC2 states x/d < 0.45 and that ratio between lengths of 
adjacent spans must not be higher than 2. Plastic method, inspired by the 
kinematic theorem of plasticity theory, sets no limits on redistribution, but 
x/d ratio must not be lower than 0.25 at any point, and high-ductility steel 
shall be used. 

Applying LR and P procedures (see above) and non-linear analysis (NL) 
to a series of two-span rectangular beams, uniformly loaded and 
characterised by the following parameters: 
width 50 cm, height 100 cm, effective height 93 cm, fy = 430 N/mm), fa = 
30 N/mm?, 

(x/d)y = 0.2; 0.225; 0.250; 0.300; 0.350; 0.400; 
1/h = 20, 25, 30, 35 
same ultimate resisting moments on support and span 
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the bearing capacities given in table 3, extended for the plastic method to 
x/d = 0.30, have been found. 


Table 3 - Two span beam (b=50 cm, h=100 cm). Parametric study for high 
values of x/d, l/h. 
Bearing capacity (kN/m). The values in brackets refer to cases with three 


hinges. 
STE GGARSER eee 
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x/d = 0.200 11.69 
|P___[36.54 


(41.34) [26.14 7A1 12.39 
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The responses in serviceability conditions as determined from non-linear 
design are described in Table 4: 
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Table 4 - Two span beam. Parametric study for high values of x/d, 1/h 
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On the results obtained, we can say what follows. 
First of all it should be pointed out that ultimate non-linear response is 


nearly always pseudo-elastic (only one plastic hinge). This is due to the type 


of dimensioning procedure adopted which delays the formation of the hinges 
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in the span. A comparison between NL and LR results is illustrated in 
graphic form in figure 2. 


(LR-NL)/NL% 





Figure 2 - Two span beam. Parametric study. Ratio (LR-NL)/NL% 


It can be seen that LR is higher than NL by nearly 10%, in very slender 
beams and for small x/d ratios. The reverse is true in opposite conditions. By 
comparing LR and P we find that, in Table 3, the two bearing capacities are 
different everywhere. The reason is that redistribution formula (1) adopted 
for LR impedes to reach the failure load simultaneously in the two critical 
zones (support, midspans). Hence, the bearing capacity P always exceeds 
LR. This occurs even when LR exceeds NL. In actual fact, as mentioned 
above, P always exceeds NL. In some cases, the difference attains 15% for 
x/d = 0.25; 20% for x/d = 0.30 (1/h = 0.35). In non-linear design, the 
phenomenon of “premature yielding” mentioned above is very common (for 
low values of x/d and 1/h). This is seen to occur in the presence of marked 
redistributions (ò close to 0.7). The bottom part of the Table also reveals the 
presence of considerable stresses in concrete. This aspect will be discussed 
later on. 

The foregoing considerations show appreciable discrepancies between 
“exact” non linear results and results provided by approximate methods. It 
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can also be noted that, in relation to serviceability conditions, non-linear 
method permits excessive redistributions. It is therefore but logical to try 
and avoid these shortcomings. 

The only available way to amend non linear EC2 method, consists of 
reducing admissible plastic rotation. Eurocode 2 admits ca 20/1000 of 
radiant up to x/d = 0.175 and decreasing values according to a hyperbolic 
law for higher values. We reach 6/1000 of radiant for x/d = 0.50. In our 
analyses we have proposed a uniform reduction of 5/1000 over the entire 
variation range. This correction improves the NL response, but does not 
fully eliminate the excess of redistribution, especially in case of pseudo- 
plastic response when the reduction of admissible plastic rotation is not very 
effective. But since non-linear design provides a full description of 
serviceability states, possible excesses in redistribution (and ensuing high 
steel stresses) can always be corrected a posteriori by a reduction of bearing 
capacity. 
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The effect of the variations made to admissible plastic rotation are shown in 
Table 5: 


Table 5 
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Table 5 shows that, with lower plastic rotation, è, values are everywhere 
increased and steel stresses maintained within acceptable range (lower than 
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0.85 fu). Yet the results are not so good in pseudo-plastic cases (see the 
beams with two hinges in table 5). 


To harmonise LR and non-linear approach with reduced plastic rotations, we 
have then proposed the following more severe redistribution formula: 


8>0.64+08®: 8>0.80 
d (2) 


which gives: è = 0.80 for x/d = 0.20, è = 1 for x/d = 0.45 


Values obtained with formula (2) are in good agreement with those of Table 
DI 


Let us however point out that the reduction of admissible plastic rotation 
increases the gap between P and LR. The difference reaches 30% for 
x/d=0.30, 1/h=25; 23% for x/d =0.25, 1/n=35. Other examples show that the 
disparity remains significant even for x/d=0.10. The only available mean to 
amend plastic results would then consist in a strong limitation of x/d. 
Practically an unacceptable restriction of its field of validity. 

We have therefore proposed to adopt for P method the same restrictions than 
for LR for what concerns the limit value of è. In actual practice to unify the 
two procedures. 


5. Problems concerning serviceability limit states 
Sa Redundant effects of cracks and imposed deformations ° 


A few remarks on an aspect of a general nature concerning crack 
formation in redundant concrete structures. 


Intuitively, it is logical to liken the effects of a crack to those of an 
imposed deformation. For instance, if in a uniformly loaded, double restraint 
beam, we introduce a hinge in the midspan section, the relative rotation 
produced by the load in the section will be equal to the rotation that would 
have to be applied by producing an imposed deformation suitable to nullify 


° F. LEVI., Cracking effects in redundant reinforced concrete structures. Acc. SC. Turin, Proceedings 
SC. Fis. 131 (1997). 
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the moment acting in the uncracked beam. It might therefore seem logical to 
attribute to redundant cracking effects the properties typical of imposed 
deformations: requirement of safety margins markedly lower than those 
assumed for loads, tendency of corresponding static conditions to undergo 
great reductions in the presence of delayed deformations. In the report cited 
above, it is demonstrated instead that likening the cracks to imposed 
deformations is legitimate only for determination of their instantaneous 
effects. In these circumstances, this approach may be useful to take 
advantage of reciprocity theorems and influence lines. It cannot be used, 
instead, to define the safety margins and in connection with delayed 
phenomena. In physical reality, the cracks are evolving imposed 
deformations which must always meet the equilibrium and compatibility 
requirements. Hence the possibility of extending to their effects the safety 
margins valid for the causes which provoked them (loads or imposed 
deformation). Same conclusion applies to influence of delayed phenomena. 


5.2 Stresses assessment at serviceability limit state 


This topic has been touched on in $ 4.1, when pointing out the close 
correlation between redistributions and steel stress. In Table 3 ($ 4.3) we 
have also briefly mentioned cases in which concrete stresses reached high 
values. Corrective measures were also suggested: reduction of plastic 
rotation, restriction of redistributions. 


Here we should like to reconsider the problem from the standpoint of a 
designer who cannot, or does not want to, perform a non-linear analysis. Our 
aim will be to supply simplified working tools for achievement of 
acceptable stress levels. 


For steel, the qualitative considerations expressed in $ 4.1 can be easily 
converted into quantitative terms if we assume a conventional ratio between 
the ultimate moment and moment at rare condition equal to 0.7 - an 
acceptable assumption for usual loading combinations. Let us calculate 
ultimate moment of a cross section attributing alternatively to materials 
design mechanical properties and characteristic properties and introducing 
the ratio X between the corresponding lever arms. We can then easily 
express the steel stress at rare level by the following simple formula: 


A) 
0, = 03 (3) 
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Assigning to Y the value 1.05, which covers with a good approximation 
the ranges within which the tension reinforcement percentage and the ratio 
between tension and compression steel sections may vary, we can write: 


0,=0.73f,, ir =0.64f,, ni 
i î (3’) 

In relation to formula (3’), the following considerations apply: 

- for determination of parameter ò,, as defined in $ 4.1, a designer always 
has on hand all necessary elements; 

- the definition of the è, parameter, characterising cracking 
redistributions, instead, would call for a non-linear analysis. It should be 
noted, however, that in rectangular beams, è; differs from 1 by less than 
10%, with maximum values for end supports of continuous beams with 
variable spans. Higher differences, of up to 15-20%, are found in T beams, 
in the same conditions. On these basis, the designer may therefore evaluate 
ò, with sufficient approximation. 

To predict concrete stresses in serviceability conditions, we may then 
follow a simple graphic procedure. 
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Let us consider fig. 3. The steeper curves, marked out with the value of 
characteristic concrete strength (fx = 20, 30,... N/mm), represent the known 
correlation between ultimate adimensional resisting moment (Ly = Muy/(bd° 
f-a) and tension reinforcement geometrical ratio (p = Ay/bd). 

The diagram is plotted for fx = 430 N/mm? 

The second family of curves (at a lesser inclination) is marked out 
instead with the o:/fwx ratio, where ©; is concrete stress in rare loading 
condition. To plot them, we have determined the serviceability moment 
corresponding to p and 0:/f-x parameters (assuming n = 15) and then divided 
this moment by 0.7 to return to the ultimate moment. In this manner, the 
second set of curves is referred to the same axes as the first set. Now, if we 
consider intersections between the curves of the first family and those of the 
second one, which correspond to a given ratio 0/fx, for instance 0.6, the 
abscissa of the intersection will show the limits of the tension steel 
percentage beyond which the stress &; = 0.6 fx is exceeded. 

Without going into further details, we should point out that the diagram 
of fig. 3 allows to take into account design redistributions (ò < 1) and 
cracking redistributions (ò, # 1). 


It is then interesting to notice that the availability of diagrams of the type 
shown in fig. 3 makes it possible to replace the notion of “admissible 
stresses” with a more realistic approach based on “admissible strength 
classes”. For example, by adopting, as in EC2, an upper limit for concrete 
stress of 0.6 fx (for exposure classes 3 - 4 - 5a - Sb), we might formulate, for 
simple bending, the following extremely simple working rules: 

- Adoption of concrete classes not lower than C 25 - 30 for tension 
reinforcement percentages of up to 1%. 

- Adoption of classes not lower than C 30 - 37 for percentages of up 
to 1.5% (the values given are conservative in order to cover the presence of 
redistributions or imposed deformations). 

It goes without saying that even “admissible strength classes” approach 
would suffer from uncertainties typically associated with the method of 
‘“admissible stresses” (for instance, choice of moduli-ratio). However, it 
avoids the reference to precise numerical limits for quantities (stresses) 
determined from purely formal calculations. 

This procedure can be extended to compression and bending. 
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6. Miscellaneous 


We shall only briefly mention a few problems discussed in the last of the 
three reports presented at the Academy, the one which is the currently being 
printed. 


6.1 Frames at the ultimate limit state 


Within the framework of the re-design of an extremely large-sized frame, 
we have compared ultimate bearing capacity determined from a linear 
analysis and from a non-linear analysis initially developed according to EC2 
method. With some surprise, the bearing capacities were found to be nearly 
the same. Indeed non-linear redistributions, localised in beams, were 
practically ineffective on columns where failure took place. In actual fact, 
the comparison was falsified by the presence of exceedingly limited 
reinforcement anchoring lengths at the beam-columns nodes. Accordingly, 
in linear approach, design process should have been halted at the loading 
level giving rise to initial slip. Working step-by-step with non-linear 
approach, instead, analysis could be carried on by imposing that relative 
displacement between reinforcement and concrete remained contained 
within ascending branch of “slip-bond” response. In these conditions, 
indeed, plastic hinges continue to work in safety conditions. 


The same structure was then analysed through the “mean values Yr” 
method (again within the framework of “slip-bond” law) and the results 
obtained were quite similar to those supplied by non-linear procedure since, 
in this case too, the results were conditioned by the poor interaction between 
beams and columns. 


6.2 Frames in serviceability conditions 


In frames, we are always in the presence of compression and bending 
stresses which cannot be handled with the procedures illustrated in fig. 3 of 
$ 5.2 for determination of concrete stresses. It proved possible, however, to 
extend the method to the presence of compression and bending with centre 
of pressure external to section by exploiting Ehler’s stratagem. This means 
that the bending moment is evaluated with reference to the tension 
reinforcement, and then the corresponding section A; is replaced by a 
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fictitious section Ay* = As - Nu/fya, where N, is the value of the ultimate 
normal force. When the centre of pressure is internal to the section, stress 
estimates require more complex graphic procedures. It is possible, at any 
rate, to make available to the designer simple elements of judgement 
enabling him to work in full safety, at the overall design and execution 
stages. 


6.3 Imposed deformations 


We shall merely summarize the conclusions of a non-linear analysis 
applied to a rather extensive series of beams and frames. 


In general, in concrete structures, at ultimate state, the incidence of 
imposed deformations that designers have normally to consider is relatively 
modest. The natural outcome is the possibility - within the framework of 
non-linear analyses - of adopting lower safety factors. For reasons 
mentioned above, the factor should even be reduced to 1 in the presence of a 
pseudo-plastic ultimate response. For pseudo-elastic responses, instead, it 
seems advisable to take into account a little model uncertainty by proposing 
a factor of the order of 1.1. Given difficulties associated with a priori 
prediction of response feature, it is logical however to adopt factor 1.1 in all 
cases, since, at all events, an increase of 10% in amplitude of imposed 
deformations has virtually no effects in pseudo-plastic conditions. 


In our third report, we also proposed an approximate design method 
founded on linear checks performed on structure assumed to be totally 
cracked (but by taking into account the variation of the tension 
reinforcement cross section along the axis). The procedure turns out to be 
fully conservative in the event of a pseudo-plastic response, and to be closer 
to reality, but still prudent, in pseudo-elastic conditions. In either case, it is 
therefore possible to adopt a safety factor as low as 1. 


In serviceability conditions, estimates are not as easy. In this case too, 
however, it is possible to obtain reliable results by working in the elastic 
field. For the application methods, the reader is referred to our reports. The 
inevitable inaccuracy of the procedure is largely compensated for by the 
intervention of delayed deformations. 
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7. Proposal of a simplified non-linear dimensioning procedure 


The method is based on a systematic use of the è, parameter, as defined 


in $4.1. 


Let us take into consideration a multi-span continuous beam subjected to 


a certain number of loading conditions. The proposal consists of performing 
the following steps: 


a): Plotting elastic diagrams of moments for a loading condition affecting 
all spans (the most logical is an evenly distributed load). 

b): Adopting the following moment redistribution criteria: 

Mr at span = Mgq on the penultimate support (for the end spans) 
M7 at span = Mg on the adjacent supports (for intermediate spans). 

c): Defining the reinforcement percentages in the critical sections and 
determining the corresponding values of òu. 

d): Re-adjusting the values of è, so as to satisfy the limits imposed by the 
redistribution formula and the need to take into account the equilibrium 
conditions 

e): Plotting elastic diagrams of moments for other loading conditions and 
adapting reinforcement percentages. 

f):  Re-adjusting values of è, as in d). 

g): Checking levels of stresses in serviceability conditions, as specified in 
$ 5.2 (Notice, however, that for steel, the preliminary checking of è, 
greatly facilitates the task). 

h): Adapting reinforcement to building requirements. 

i): If deemed necessary, checking cracking and deformation conditions by 
means of simplified methods. 
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Time-Optimal Climb of Next-Generation Jet Fighter 
Aircraft: Classical and Computational Results! 


A. MIELE” and J. B. DABNEY 


Abstract. This paper describes an investigation of the climb 
performance of a jet fighter aircraft capable of maneuvering at angles 
of attack in the post-stall region. The flight dynamics model accounts 
for all of the factors necessary to accurately characterize aircraft 
motion. The maneuvers studied consist of climbs from 5 to 35 kft and 
include two classes: climbs in a purely vertical plane and climbs in 
which the bank angle is allowed to change. 

The main conclusion is that time-optimal climb maneuvers do not 
require flight in the post-stall region. Instead, time-optimal climbs 
require relatively low angles of attack and are characterized by three 
distinct segments: a central segment often flown with a load factor of 
nearly 1 and two terminal segments (dive or zoom) to and from the 
central segment. The central and final segments are nearly 
independent of the initial conditions, instead being dominated by the 
final conditions. In all cases, the initial segment of the maneuver 
consists of a zero load factor dive so as to reach the central segment. If 
the bank angle is allowed to change, certain trajectories with excess 
speed at the end of the zoom maneuver conclude with a moderately 
high load factor maneuver to align the velocity vector with the horizon. 

Comparison of two-dimensional and three-dimensional trajectories 
shows that, in the case of low final speed, there is a modest benefit in 
flight time (about 5 percent) by performing a 180 deg roll prior to 
reaching the desired final conditions. 
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1. Introduction 


Time-optimal climbs are important maneuvers for jet fighter aircraft. For 
example, the interception of a high-altitude aircraft requires the interceptor 
to climb to an appropriate altitude as rapidly as practical, and to arrive at 
that altitude in a state consistent with the capability to maintain level flight. 
As another example, it may be desirable for an aircraft, after expending 
ordnance against a ground target, to rapidly climb above enemy defenses, 
and to arrive at a higher altitude with a speed such that the ability to 
maneuver evasively is preserved. As a consequence of the importance of 
time-optimal elimbs of jet fighter aircraft, several studies of these 
trajectories have been reported over the last fifty years (Refs. 1 - 11). 

This paper describes a numerical study of time-optimal climbs for next- 
generation jet fighter aircraft capable of maneuvering at high angles of 
attack. The results of the computational study are compared with those 
produced via the classical theory of optimal climb formulated in Refs. 2 - 3. 

The optimal trajectories computed here serve several purposes. First, 
they provide insight into the actual usefulness of the ability to maneuver at 
extreme angles of attack. Second, they serve as basis for guidance 
trajectories that approximate the optimal climb. Third, they illustrate the 
advantages accruing in certain’ situations via three-dimensional 
maneuvering, rather than via the purely two-dimensional maneuvering 
studied in Refs. 1 - 11. 

The hypothetical aircraft used for this study is described in detail in Ref. 
12. Briefly, the model characterizes a hypothetical jet fighter aircraft for 
flight speeds ranging from low subsonic (Mach number less than 0.1) to 
moderate supersonic (Mach number approaching 2.0), altitudes from sea 
level to 82 kft, angles of attack from 0 to 90 deg, and parametrized wing 
loading and thrust-to-weight ratios. 


2. System Description 


In this section, we describe briefly the mathematical model of a jet 
fighter aircraft flying in a vertical plane. For a more detailed discussion of 
the aircraft model, refer to Ref. 12. We make the following assumptions: 
(A1) flight takes place at altitudes from sea level to 82 ft; (A2) flight takes 
place over a flat Earth; (A3) there is no wind; (A4) weight change during 
climb is neglected; (A5) aircraft is controlled via angle of attack: (A6) 
engine is operated at maximum thrust; (A7) thrust is directed along the 
aircraft reference line (thrust vectoring is not considered); hence, thrust 
angle of attack is the same as aerodynamic angle of attack. 
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2.1. Differential System. With the above assumptions, the motion of the 
aircraft is described by the following differential system for the along-track 
distance x, altitude h, velocity V, and flight path angle Y: 


x= Vcosy, (la) 
y= Vsiny, (1b) 


V=(Tcosa - D- Wsiny)/m, 
(le) 

y=(L+ Tsina — Wcosy)/mV, (1d) 
in which the dot denotes derivative with respect to the time t. Here, 0<t £ 1, 


where t is the final time. The quantities appearing on the right-hand side of 
(1) are the thrust T, drag D, lift L, weight W, mass m, and angle of attack a. 


2.2. Functional Relations. In the system (1), the weight W is assumed to 
be constant, 


W=mg,, (2) 
which is the product of the constant mass m and the sea-level acceleration of 


gravity g.. Generally speaking, the aerodynamic forces can be represented 
by the functional relations 


D=D(h,V,0), (3a) 

Le 1lh;V0). (3b) 
For a given power setting, the engine thrust is described by the functional 
relation 


T=T(h;V). (4) 


2.3. Inequality Constraints. Inspection of the system (1) in light of (2) - 
(4) shows that the time history of the state x(t), h(t), V(t), y(t) can be 
computed by forward integration for given initial conditions, given control 
o(t), and given final time ©. In turn, the control is subject to the two-sided 
inequality constraint 


O0<o<7T/2 (5a) 


or 
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0<Sas<a,, (Sb) 


which must be satisfied everywhere along the interval of integration. 
Inequality (Sa) applies if the aircraft is permitted to operate in the 
supermaneuver range; Ineq. (Sb) applies if the angle of attack is forced to be 
less than that for maximum lift coefficient o». Regardless of whether (5a) or 
(Sb) applies, this inequality constraint can be accounted for via a 
trigonometric transformation discussed in detail in Ref, 12. 


2.4 Atmospheric Model. The atmospheric density p(h), temperature 
8(h), and speed of sound a(h) are computed using approximations to the 
1962 U.S. Standard Atmosphere (Ref.13) as described in Ref. 14. Details of 
the atmospheric model are presented in Ref. 12. 


2.5. Aerodynamic Forces. The aerodynamic forces (3) are given by 
D'<%/2)CgpSV°, (6a) 
L'=(1/2)C psv”, (6b) 


where Cp is the drag coefficient, C, the lift coefficient, and S a reference 
surface area. The aerodynamic coefficients represent a hypothetical fighter 
aircraft based on the data of Refs. 14 - 16. The coefficients used here 
account for angles of attack up to 90 deg, and also take into account the 
dependence on Mach number of the maximum lift coefficient. 

In the flight regimes of interest, the coefficients of drag and lift can 
be represented by the functional relations (Ref. 12) 


Cp =Cp(0,M), (7a) 
C, =C,(0,M), (7b) 
where M is the Mach number, with 


M=V/a(h). (8) 
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2.6. Engine Model. The engine thrust approximation is split into two 
regions: troposphere and stratosphere. The engine thrust (4) is approximated 
as 


T=T_,F,(h,M), (0) 


ref 
where T,x is the maximum sea-level static thrust and Fy(h,M) is a 
normalized thrust function. In the troposphere, Fr is computed using 


polynomial approximations to the ratio T/T,x for the J85 turbojet engine 
(Refs. 12 and 17). In the stratosphere, Fy is represented as 


F, (h,M)=(p/p.)F;(h.,M), (10) 


where p- is the density at the tropopause altitude h+. For details, see Ref. 12. 


2.7. Path Constraint. To avoid structural damage to the aircraft, the 
aerodynamic lift is constrained according to the relation 


L+Tsino<n, W (11) 


max ? 


where nmax is the load factor limit. This constraint can be imposed via a 
slack variable transformation (Refs. 18 - 19). Due to the high load factor 
limit of the hypothetical jet fighter aircraft [see (12c)], Ineq. (11) is not 
violated in any of the computer runs performed for this study. Therefore, 
Ineq. (11) can be omitted from the computations. 


2.8. Boundary Conditions. The initial conditions (t = 0, subscript i) are 
inputs to the optimization problem. The final conditions (t =t, subscript f) 
are also inputs to the optimization problem. 


2.9. Supplementary Data. The following data have been used in the 
numerical experiments: 


reference surface area S = 237.0 ft°, (12a) 
aircraft nominal weight W = 10 kIb, (12b) 
load factor limit nmax = 9, (12c) 


maximum engine thrust with afterburning T;e:=12.0 KIb, (12d) 
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maximum engine thrust without afterburning T;e = 6.8 KIb, 

(12e) 

sea-level acceleration of gravity g. = 32.174 ft/sec? . (12f) 
Clearly, the wing loading is relatively low, W/S = 42.2 Ib/ft.° On the other 
hand, the thrust-to-weight ratio is relatively high, Tref/W = 1.20 with 
afterburning and Tref/W = 0.68 without afterburning. 


3. Optimization Problem 


The minimum-time climb problem (Problem P1) can be formulated as 
follows: 


minI=7, (13) 


where t is the final time. The unknowns include the state variables x(t), h(t), 
V(t), y(6) the control variable o(t), and the parameter t. 


3.1. Time Normalization. In treating problem (13), it is convenient to 
normalize the dimensional time t via the transformation 


n=t/T (14a) 


where n is the dimensionless time and 1 the final time. Clearly, 
VENITE SET (14b) 


3.2. Approach. The optimization problem under consideration is of the 
Mayer type and can be solved via the sequential gradient-restoration 
algorithm for optimal control problems (SGRA, Refs. 18 - 19). The 
approach taken is to produce manually a nominal control history that results 
in a solution of the differential system (1), satisfying the specified initial 
conditions and the path constraint (11), and approximately satisfying the 
final conditions. This solution is used as the starting nominal trajectory for 
SGRA; this algorithm produces a sequence of feasible trajectories, each 
characterized by a lower value of the performance index I than the previous 
feasible trajectory. The algorithm terminates when the constraint error and 
the error in the optimality conditions become smaller than certain 
preselected small positive tolerances. 


3.3. Computer Runs. For the two-dimensional case, 6 optimum climb 
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trajectories (Trajectories TI - T6) were computed by considering maximum 
engine thrust with afterburning [see (12d)] and various combinations of the 
following parameter values: 


x(0)= O kft, x(t) =ifree, (15a) 

h(0)= Skft, ha) (15b) 

V(0)= 0.5kft/sec, V(t)=05, 1.0, 1.5 kft/sec, (15c) 

Y(0)=0, Y()=0, (15d) 
or 

Y(0) = free, Y(t) = free. (15e) 


State variables which are not prescribed must be considered as free and must 
be determined from the optimization process together with the main 
unknown, the final time. The computations via SGRA were done with an 
Intel Pentium II processor using the C-language and double-precision 
arithmetic. 

For convenience, Table 1 presents a summary of the trajectory 
parameters identifying Trajectories T1 - T6 [see (15)]. Trajectories T1 - T6 
are plotted in Figs. 1 - 6, along with a tabulation of the key data and results. 


Table 1- Two-dimensional optimal climb trajectories. 


Trajectory V(0) [kft/sec] V(1) [Kft/sec] Y(0) [deg] Y(© [deg] t[sec] 


[1 0.5 0.5 0 0 81.1 
Te 0.5 1.0 0 0 90.8 
23 0.5 1.5 0 0 159 
T4 0.5 0.5 DI 47.5 74.9 
T5 0.5 1.0 -18.1 DE 90.2 
T6 0.5 Hei -25.8 2.4 113.3 





4. Classical Theory of Optimal Climb 


It is of interest to compare the trajectories produced in this study with the 
classical theory of optimal climb of jet aircraft flying in a vertical plane 
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(Refs. 1 - 5). As developed in Refs. 2 - 3, the classical theory is based on the 
following assumptions: (i) the thrust is nearly tangent to the flight path; (ii) 
the normal acceleration is neglected by comparison with the tangential 
acceleration; and (iii) the drag is computed approximately via the simplified 
assumption that the load factor is 1, hence L = W. Because hypothesis (iii) 
implies 0 = o.(h, V), we can rewrite relations (3a) and (4) in the form 

T=T(h,V), D=D(h,V). (16) 
Therefore, the excess power per unit weight, 

E=(T-D)V/W, (17a) 


is also a function of the form 


E=E(h,V); (17b) 


for that matter, the so-called omega function (Refs. 2 - 3), 


O=—--——_, (18a) 


has also the form 
© = ®(h, V). (18b) 
Note that the E-function has the dimension of a velocity, while the @- 
function is dimensionless. In light of the above hypotheses, the flight path 
angle y replaces the angle of attack © as the control and the minimum-time 
trajectory consists of the following segments (Refs. 2 - 3): 
y=-1/2, (19b) 


Y=+7/2. (19c) 
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Thus, the classical theory suggests that the optimal climb trajectory consists 
of three segments: Equation (19a) describes the central segment of the 
optimal climb trajectory; Equations (19b) and (19c) describe the initial and 
final segments. Depending on the boundary conditions, the initial segment 
can be a vertical dive or zoom; analogously, the final segment can be a 
vertical dive or zoom. With reference to the (V, h)-domain, the choice of 
(19b) vs. (19c) for the initial and final segments of the optimal trajectory 
depends on the location of the initial and final points with respect to the line 
0=0. 

For the purpose of comparison, the trajectory obtained via the classical 
theory of optimal climb (Refs. 2 - 3) is plotted in the (V, h)-domain (dashed 
line, Figs. 1 - 6) together with the optimal trajectory computed via SGRA 
(solid line, Figs. 1 - 6). 


4.1. Remark on Altitude. Assume that there is an altitude constraint 
having one of the three forms below, 


he hL. (20a) 
hsho (20b) 
h, <h<hy, (200) 


namely, the aircraft is required to fly in a bounded altitude domain; indeed, 
the subscripts L and U stand for lower bound and upper bound, respectively. 
Then, the optimal trajectory might include, in addition to (19), the subares 

ha hi Ò (21a) 
and/or 

h=h, (21b) 
depending on the boundary conditions of the problem. Subarcs of type (21) 
were found in Refs. 20 - 21; they occur in this work because of the altitude 
constraint 


h20. (22) 


See the dashed trajectories in Figs. 1 - 6. 
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5. Analysis 


In this section, we discuss the computed optimal trajectories in light of 
the classical theory of optimal climb trajectories for turbojet aircraft (Refs. 2 
- 3). Inspecting Trajectories T1 - T6, it is evident that, as predicted by the 
classical theory of optimal climb, each trajectory consists of three segments. 
The central segment is a relatively long duration climb. The initial segment 
consists of transferring the aircraft from the initial conditions to the central 
segment; the final segment is a low-angle-of-attack zoom maneuver leading 
from the central segment to the specified final conditions. 


5.1. Central Segment of Optimal Trajectory. Examining Figures 1 - 6, 
it is evident that the optimal climb trajectories involve a central segment 
flown with a load factor of approximately 1, thereby justifying the 
assumption L = W employed in Refs. 2 - 3. In this section, we discuss the 
influence of the final speed and terminal flight path angles on the central 
segment. With reference to the (V,h)-domain, we also compare the 
computed central segment with the trajectory © = 0 [Eq. (19a)] predicted by 
the classical theory of optimal climb (Refs. 2 - 3). 

Figures 1 - 6 illustrate the computed optimal climb trajectories for three 
final speeds [see (15c)] and for terminal flight path angles either zero [see 
(15d)] or free [see (15e)]. Each figure also shows the corresponding 
classical trajectory predicted via (19). Examining these trajectories, it is 
clear that the influence of the final speed on the central segment is quite 
strong. 

To compare the computed optimization results to the classical theory, 
consider again Figs. 1 - 6. The central segment of the trajectories featuring 
higher final speeds (1.0 and 1.5 kft/sec) follows relatively closely the © = 0 
trajectory. On the other hand, the central segment of the trajectories 
featuring a lower final speed (0.5 kft/sec) bear less resemblance to the © = 0 
trajectory. 


5.2. Initial Segment of Optimal Trajectory. All of the trajectories 
presented here begin at a low initial speed. In each case, the initial segment 
of the optimal trajectory consists of an acceleration at low angle of attack 
followed by a pitch-up maneuver with low load factor to intercept the 
central segment of the optimal trajectory. 

Examining Figs. 1 - 6, it is evident that relatively steep dives are 
employed, albeit not the vertical dives predicted by the classical theory. The 
initial acceleration maneuver results in a dive angle of approximately 15 - 30 
deg. Examining Table 1, it is evident that, even for initial flight path angle 
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free, the optimal trajectories do not involve vertical dives. 


5.3. Final Segment of the Optimal Trajectory. The final segment of 
the optimal trajectories consists of a zoom maneuver at low load factor so as 
to reach the specified final boundary conditions. The vertical zooms 
predicted by the classical theory of optimal climb are not evident although 
relatively steep climbs are employed, particularly in the case of low final 
speed. The magnitude of the benefit accrued by allowing the final flight path 
angle to be free depends on the desired final speed. In the case of low final 
speed (Trajectory T4 of Table 1), the optimal final flight path angle is Y(1) = 
47.5 deg and there is a slight time benefit (approximately 8 percent) with 
respect to the case where y(t) = 0 (Trajectory T1 of Table 1). On the other 
hand, in the cases involving higher final speed, the optimal final flight path 
angle Y(t) is small positive, and there is no benefit with respect to the case 
where Y(T) = 0. 


6. Three-Dimensional Maneuvering 


The next generation of jet fighter aircraft can be expected to have the 
capability to generate very high rates of roll. To assess the importance of 
this capability within the frame of optimal climb maneuvers, three additional 
optimal trajectories were computed using a three-dimensional model and 
hence equations of motion of higher order. 
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Figure IA 
Optimal climb trajectory Tl: low final speed, vertical plane, fixed terminal flight 
path angles. 


h(0) = 5 kft, V(0) = 0.5 kft/sec, Y(0) = 0.0 deg. 
h(1) = 35 kft, V(1) = 0.5 kft/sec, YT) = 0.0 deg. 
t=81.l'sec. 





Figure 1B- Optimal climb trajectory T1: data and results. 
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Figure 2A 
Optimal climb trajectory T2: medium final speed, vertical plane, fixed terminal flight 
path angles. 





h(0) = 5 kft, V(0) = 0.5 kft/sec, Y(0) = 0.0 deg. 
h(1) = 35 kft, V(t) = 1.0 kft/sec, YT) = 0.0 deg. 
T=90.8 sec. 





Figure 2B- Optimal climb trajectory T2: data and results. 
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Optimal climb trajectory T3: high final speed, vertical plane, fixed terminal flight 


path angles. 


h(0) = 5 kft, V(0) = 0.5 kft/sec, 
h(q) = 35 kft, V(1) = 1.5 kft/sec, 
T= 113.9 sec. 


Y(0) = 0.0 deg. 
Y(T) = 0.0 deg. 





Figure 3B- Optimal climb trajectory T3: data and results. 
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Optimal climb trajectory T4: low final speed, vertical plane, free terminal flight path 
angles. 


h(0) = 5 kft, V(0) = 0.5 kft/sec, Y(0) = 9.3 deg. 
h(a)= 35 kft, V(1) = 0.5 kft/sec, Y(T) = 47.5 deg. 
T= 74.9 sec. 





Figure 4B- Optimal climb trajectory T4: data and results. 
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Optimal climb trajectory T5: medium final speed, vertical plane, free terminal flight 
path angles. 


h(0)=5.kft, V(0) = 0.5 kft/sec, Y(0) = -18.1 deg. 
h(t) = 35 kft, V(1) = 1.0 kft/sec, Y(T) = 5.2 deg. 
T=90.2 sec. 





Figure 5B- Optimal climb trajectory TS: data and results. 
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Optimal climb trajectory T6: high final speed, vertical plane, free terminal flight path 
angles. 


h(0) = 5 kft, V(0) = 0.5 kft/sec, Y(0) = -25.8 deg. 
h(1) = 35 kft, V(1) = 1.5 kft/sec, Y(T) = 2.4 deg. 
ti=113.3'scc: 





Figure 6B- Optimal climb trajectory T6: data and results. 
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6.1. Differential System. To study maneuvering in three dimensions, we 
replace the fourth-order differential system (1) with the sixth-order 
differential system (23) by including the cross-track distance y and velocity 
heading angle Y as state variables, 


X= Vcosycosy, (23a) 
y= Vceosysiny, (23b) 
h= Vsiny, (23c) 
V=(Tcosa —- D- Wsiny)/m, (23d) 
y=|[(L + Tsina.)cosy — Wcosy]/mV, (23e) 
X=(L+Tsino.)sinu / mVcosy. (23f) 


Here, the bank angle u(t) is a control in addition to the angle of attack ot). 


6.2. Computer Runs. Three additional computer runs were performed 
using this three-degree-of-freedom aircraft model. The resulting Trajectories 
T7 - T9 use the same parameter combinations as Trajectories TI - T3. For 
convenience, Table 2 presents a summary of the trajectory parameters 
identifying Trajectories T7 - T9 [see (15)]. Trajectories T7 - T9 are plotted 
in Figs. 7 - 9, along with a tabulation of the key data and results. 


Table 2- Three-dimensional optimal elimb trajectories. 

Trajectory V(0) [kft/sec] V(7) [kft/sec]  y(0) [degl y(© [deg] t[sec] 
s 0) 0.5 (05 0 0 77.6 
T8 0.5 1.0 0 0 90.8 
T9 0.5 15 0 0 113.9 


6.3. Analysis. Examining Trajectories T7 - T9, it is evident that three- 
dimensional optimal climb trajectories also consist of three segments. The 
initial and central segments are identical to the initial and central segments 
of the corresponding two-dimensional optimal climb trajectories. The final 
segment begins with a rapid roll to a bank angle of 180 deg. In the case of 


TIME-OPTIMAL CLIMB OF NEXT-GENERATION JET FIGHTER AIRCRAFT 165 


low final speed, the final segment ends with a turn in a near-vertical plane at 
relatively high angle of attack so as to decrease rapidly the flight path angle 
and achieve the specified final conditions. In the cases involving higher final 
speeds, the final segment ends with a low-angle-of-attack zoom to reach the 
specified final conditions. 

In the case of low final speed (Trajectory T7), relative to the 
corresponding vertical plane trajectory (Trajectory T1), there is a modest 
benefit in flight time (5 percent reduction) accrued because of the final 
rolling maneuver. For the cases involving higher final speeds (Trajectories 
T8 and T9), relative to the corresponding vertical plane trajectories 
(Trajectories T2 and T3), there is no benefit accrued from the final rolling 
maneuver; indeed, in both the two-dimensional and the three-dimensional 
cases, the final segment is flown at very low angle of attack. 
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Optimal climb trajectory T7: low final speed, roll permitted, fixed terminal flight 
path angles. 


h(0) = 5 kft, V(0) = 0.5 kft/sec, Y(0) = 0.0 deg. 
h(t) = 35 kft, V(1) = 0.5 kft/sec, YT) = 0.0 deg. 
T=77.6 sec. 





Figure 7B- Optimal climb trajectory T7: data and results. 
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Optimal climb trajectory T8: medium final speed, roll permitted, fixed terminal flight 
path angles. 


h(0) = 5 kft, V(0) = 0.5 kft/sec, Y(0) = 0.0 deg. 
h(t) = 35 kft, V(1) = 1.0 kft/sec, YT) = 0.0 deg. 
T=90.8 sec. 
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Figure 8B- Optimal climb trajectory T8: data and results. 
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Optimal climb trajectory T9: high final speed, roll permitted, fixed terminal flight 
path angles. 


h(0) = 5 kft, V(0) = 0.5 kft/sec, Y0) = 0.0 deg. 
h(t) = 35 kft, V(1) = 1.5 kft/sec, YT) = 0.0 deg. 
t= 113.9 sec. 
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Figure 9B- Optimal climb trajectory T9: data and results. 
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7. Remark on Afterburning 


The minimum-time trajectories discussed in this paper were computed 
using maximum thrust with afterburning. Sometimes, it might be desirable 
to climb with afterburner off if the intent is to contain fuel consumption 
instead of time. 

As an example, consider a climb from h(0) = 5 kft to h(t) = 30 Kft, with 
V(0) = V(1) = 0.5 Kft/sec and terminal flight path angles Y(0) and yY(T) free. 
Reducing the reference thrust from Ter = 12.0 kIb (afterburner on) to Ter = 
6.8 KIb (afterbumer off) is accompanied by the following effects: (a) flight 
time becomes nearly three times larger; (b) fuel consumption per unit time 
becomes nearly four times smaller; (c) fuel consumption during climb 
decreases by approximately 20 percent. 


8. Conclusions 


Generally speaking, optimal climb trajectories for a jet fighter aircraft 
capable of maneuvering at high angles of attack do not require flight at high 
angles of attack. The optimal climbs require relatively low angles of attack 
and are characterized by three distinct segments: a central segment often 
flown with a load factor of nearly 1 and two terminal segments (dive or 
zoom) to and from the central segment. lf roll is permitted, certain 
trajectories with excess speed at the end of the zoom maneuver conclude 
with a moderately high load factor maneuver to align the velocity vector 
with the horizon. In all cases, the initial segment of the maneuver consists of 
a zero load factor dive followed by a pitch-up so as to reach the central 
segment. Comparison of two-dimensional and three-dimensional trajectories 
shows that, in the case of low final speed, there is a modest benefit in flight 
time (about 5 percent) by performing a 180 deg roll prior to reaching the 
desired final condition. 

The trajectories presented here can serve at the basis for a subsequent 
study to develop guidance trajectorics that will approximate the optimal 
trajectories. Further trajectory optimization studies are needed to 
characterize optimal climbs that include changes in heading. It is also 
desirable to study climbs in which time-to-intercept a high speed intruder is 
minimized. 
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Fracture mechanics and fractal morphologies: 
relevance of the non-integer scale-exponents 


Alberto CARPINTERI" 


Abstract. Fractal concepts are used to study the complex shapes of 
fracture surfaces, as well as damage phenomena presenting statistical 
characteristics. While at the beginning of the loading process the 
microcracks can be considered as two-dimensional surfaces, as the 
load is increased the microcracks grow, coalesce and form an invasive 
fractal set with a dimension larger than two. When the evolving 
dimension assumes the notable value of 2.5, the microcrack set may be 
considered as extremely disordered, and percolation of the cracks is 
very likely, in particular for sufficiently large specimens. Percolation 
favours catastrophic and brittle behaviours and produces fracture 
surfaces of a Brownian character, where the local dimension is again 
2.5. An additional assumption is that of considering material resisting 
sections at peak stress as lacunar fractals of a local dimension equal to 
1.5. The above-mentioned fractalities produce the well-kKnown scale 
dependence of fracture energy and tensile strength, respectively. 


Riassunto. / concetti della geometria frattale vengono utilizzati per 
analizzare le morfologie complesse delle superfici di frattura, così 
come i fenomeni di danneggiamento che presentano caratteri statistici. 
Mentre all'inizio del processo di sollecitazione le microfessure possono 
essere considerate come superfici bidimensionali, al crescere della 
sollecitazione le microfessure crescono, coalescono e formano un 
dominio frattale invasivo con dimensione maggiore di due. Allorché la 
dimensione nella sua evoluzione assuma il valore notevole 2,5, 
l'insieme delle microfessure può considerarsi estremamente 
disordinato, così che la frattura completa diventa probabile, in 
particolare nel caso di elementi sufficientemente grandi. La frattura 
completa coinvolge comportamenti meccanici catastrofici, e quindi 
fragili, e produce superfici libere di tipo Browniano, ove la dimensione 
locale si rivela essere nuovamente pari a 2,5. 

Una ulteriore ipotesi è quella che considera la sezione resistente del 
materiale alla sollecitazione di picco come un dominio frattale 
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lacunare con dimensione locale pari a 1,5. Le suddette caratteristiche 
frattali producono le ben note dipendenze dalla scala, rispettivamente 
dell'energia di frattura e della resistenza a trazione. 


1. Introduction 


As a new mathematical tool, Fractal Geometry [1-3] has recently been 
used to describe irregular phenomena in various scientific fields. In addition 
to their original use for measuring the length of irregular coastlines, it has 
been shown that, in the area of the strength and fracture behaviour of 
engineering materials, fractals are able to represent the fracture surfaces in 
two- or three-dimensional problems [4, 5]. Up to now, almost all research 
work on the fractal aspects of fracture mechanics of disordered materials has 
been concentrated on the evaluations of the fractal dimensions of fracture 
trajectories or surfaces. On the other hand, fractal geometry is not only 
useful for studying complex shapes, but may also prove a powerful tool in 
describing damage phenomena with statistical characteristics [6, 7]. For 
example, in disordered materials like rocks, concrete and ceramics, there are 
many random parameters (i.e., position, size and orientation of the pre- 
existing microcracks). The microcracks propagate during the loading 
process, and at the same time the parameters evolve according to the 
statistical character of the microcracks. 

In this paper, fractal concepts will be used to evaluate the evolving 
degree of disorder in the crack size distribution of the material. A fractal 
dimension measuring the microcrack size distribution is defined, and the 
variation of this parameter is continuously computed during the microcrack 
propagation process. It is found that the defined fractal dimension increases 
with the propagation of the cracks. In other words, it turns out that the 
materials become progressively more disordered with the development of 
the crack net, from the initial loading up to the failure stage. 

Whereas at the beginning of the loading process the microcracks can be 
considered isolated and the topological dimension of the set is two, i.e., a 
discontinuous surface in three-dimensional space, as the load is increased 
the microcracks grow, coalesce and form a more complex set. This set may 
be described and represented as an invasive fractal set with a dimension 
larger than two and smaller than three. When the evolution of the rupture 
phenomenon produces a dimension beyond the threshold value 2.5, the 
microcrack set may be considered as extremely disordered [8-10], and 
global failure, which is represented by the percolation of the cracks, 
becomes very likely. In this case, in fact, the size of the largest crack is more 
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than proportional to the specimen size, and for sufficiently large specimens 
percolation prevails. In addition, we can state that even the final fracture 
surface, after percolation and separation of the specimen into two pieces, 
presents a local fractal dimension equal to 2.5 [8-10]. In other words, it is a 
Brownian surface. In this connection, it is interesting to consider that, as 
large specimen sizes favour microcrack percolation, so at the same time they 
present more catastrophic and therefore brittle mechanical behaviours [11- 
15l: 

An additional assumption may be that of considering material ligaments 
(at peak stress) as lacunar fractals of a dimension comprised between 1 and 
2 [8-10]. If the experimental investigation is carried out over more than one 
order of magnitude in the scaling range, it is clear how not just a single and 
unique fractal dimension of fracture surface or material ligament may be 
measured by varying the structural scale. This apparent anomaly is in fact 
due to the different fractal character which emerges at different scales. In 
nature there are no ideal geometrical fractals; natural morphologies are 
generally random multifractals due to the limited size of the heterogeneities 
or perturbations (e.g., concrete aggregates). Whereas ideal fractals do not 
present any characteristic length, for random multifractals it is possible to 
distinguish between small scales, at which disorder prevails, and large 
scales, at which order prevails. For fracture surfaces, the microscopic 
disorder is that of Brownian surfaces with fractal dimension equal to 2.5, 
whereas the macroscopic order is that of Euclidean surfaces, with integer 
dimension equal to 2. On the other hand, for material ligaments at peak 
stress, the microscopic disorder is that of a lacunar surface with fractal 
dimension equal to 1.5, whereas the macroscopic order is that of a classical 
cross-sectional plane with integer dimension equal to 2 [16-21]. 

In the framework of this global fractal view, when the scale of 
observation tends to zero, the nominal tensile strength tends to infinity, 
whereas the nominal fracture energy tends to zero. If the purpose is that of 
defining constant material properties (strength or toughness), also at the 
microscopic scales, renormalization group theory allows us to find 
anomalous quantities with physical dimension equal to that of the stress- 


- È È -3R . ; 
intensity factor, 1.e., [F][L] . It emerges that, at the microscopic scale, 


strength and toughness lose their macroscopic identity and converge to the 
same intermediate quantity [8-10]. 
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2. Microcrack Coalescence and Percolation 


2.1 Defect Size Distribution of Self-Similarity 


In the present section, two-dimensional structures with a multitude of 
cracks or voids of a given size-distribution are considered. Three hypotheses 
are assumed [6, 7]: 

(1) the structure is macroscopically homogeneous; 
(2) the structure is macroscopically isotropic; 
(3) the interaction between the imperfections is negligible. 


As a first case, let us consider a set of similar structures, where a 
multitude of cracks and/or polygonal voids of constant size a are embedded 
(Fig. 1.a). They can be considered as specimens of the same material, and 
their failure occurs when the fracture condition at the imperfection of the 
most critical orientation is reached. Since such imperfections are all of the 
same size, we can state that no size effect is present in this ideal case. 

Let us now consider a set of similar structures, where the imperfections 
have a constant size which is proportional to the size of the structure (Fig. 
I.b). In this case, they cannot be considered as specimens of the same 
material. The failure occurs at the imperfection of the most critical 
orientation which is of a size proportional to the size of the structure. Since 
the influence of the other imperfections is assumed negligible, this case is 
altogether equivalent to that of a unique imperfection of size proportional to 
that of the structure. Thus, the strength size effect is represented by a linear 
In versus In b diagram with slope —@, where o. is the power of the stress 


singularity produced by the imperfection. If we have imperfections of the 
same size, but with different shapes (e.g., cracks and circular pores), the 
fracture condition has to be evaluated only for the imperfections with the 
maximum value of 0 (0<0.< 12). 

Let us consider finally a set of similar structures where the imperfections 
with the most dangerous shape (a = lic) have a probability density p(a) 
of size distribution (Figs. l.c and 2.a). We can state that if the size 
distribution p(a) is such that the maximum size a,,,x is proportional to the 
linear scale b, then the strength size effect will be represented by a linear In 
G; versus In b diagram with slope —©&,mx. The above hypothesis is very 


restrictive and is valid only when the probability density of size distribution 
p(a) presents particular properties. If p is the density of the imperfections 
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(number of imperfections per unit area), the maximum size a, can be 
defined as follows [6, 7]: 


Pb? Pd rm) = ARR IE 1) 


where the factor (1/27) is due to the fact that all the imperfection orientation 
angles © are equally likely. If a geometrically similar structure of 
characteristic size kb is considered and the above hypothesis is assumed 
(Fig. 1.c), we can write 


1 
P(Kb) Pam) Dr Aa Am = 1. (2) 


Since a, iS a function of p and b, as well as of the ranges Aa and Ao, it 
follows that Egs. (1) and (2) must be valid for any defect size a 


pb°p(a) Ri Aa AQ = 1, (3.4) 


pk°b°p(ka) E Aa AQ = 1. (3.b) 
From Egs. (3.a) and (3.b) it follows that 
p(a)=k°p{ka), Va>>a, VkeR', (4) 
and then function p assumes the form 
e 
Pal= sz» Va>>a, (5) 


where C is a constant with the physical dimension of a length, and a is the 
average defect size. Equation (5) will be referred to as the defect size 
distribution of self-similarity. The related cumulative distribution function P 
is 
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Pa) = f(x) &x= ("° px) + [" È ax, (6) 


do X 


where 4, >>a is the value beyond which the descending branch of the 


function p can be approximated by Eq. (5). Performing the integrations in 
Eq. (6), we obtain 


CI Ù ve 
Pla) = Po +|- | =P+t_—--. (7) 
Xda Ugo |;@ 
Since for d + ©© the cumulative distribution P(a) > 1, 
i 
P+t_-=1, 
fue ang (8) 
and then (Fig. 2.b) 
1% 
P(a)=1--, for a>a, (9) 
a 
with 


Generally speaking, when the condition of self-similarity is not satisfied, 
the cumulative distribution function P will have the following form: 


P(a)= 1- E VIRESrrUpsirRa. (11) 
a 
with 
Caldo (12) 


In this case, the strength size effect can be represented by a linear In 6; 
versus In b diagram with slope —A, with 
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o(Y) 


My, N)= st 


’ (13) 


where o is the power of the stress-singularity produced by a re-entrant 
corner of amplitude Y, and the exponent © depends on the secondary 


features of the material (e.g., density of imperfections, size distributions of 
the less dangerous defects, etc.). The probability density of size distribution 
in the general case — see Eq. (11) — is 


dP C 
Pa) = = ol for a> 49, (14) 


which becomes Eq. (5) when N = 1. 
The power-law Ramberg-Osgood material with the hardening stress- 


strain relation 
n 
€ (0) 
— [= Nli<xait<ic07 (15) 
€0 So 


may also be applied to study the state of affairs near a re-entrant corner with 
angle Y. Hutchinson [22] and Rice and Rosengren [23] gave the stress- 


singularity at the crack tip for this material as 
1 


= : (16) 
n+1 


with % as the exponent of r, i.e., r*. When n = 1, the linear elastic 
behaviour A = 1/2 is refound. When n> e and A+ 0, we have a rigid- 


perfectly plastic material. The stress-singularity then vanishes and the 
plastic stress-intensity factor attains the physical dimensions of stress [6]. 
This limit case is analogous to a re-entrant angle with y+> in linear 


elasticity [24]. For y > 0 and n > I, the symmetrical stress field at the re- 
entrant corner tip may be described as 


GERI O), (17) 
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in which the power À is given by 


1 


My.n)=2o(Y) nh 


(18) 


If the single-crack analysis given above is extended to the case of many 
defects of a given size distribution in nonlinear materials, the strength size 
effect can be represented by the Ino, versus In b line with a slope of -Z. 
From Egs. (13) and (18), it is found that 


XY, n, Nj=t— ; 
(Y ) w (19) 
The two-dimensional treatment can be extended to three-dimensional 


structures with polyhedral cavities. Equation (19) can be generalized as 


_ 20(9) 
MO =) (20) 


in which (N — 1) instead of N has been written, because Egs. (3) in three 
dimensions yield 


p(a)=k°p(ka), Va>>a, VkeR', (21) 
and hence 
(& 
P(a)=-, (22.a) 
a 
G 


(22.b) 


prevail for sufficiently large a. Here, © represents a re-entrant solid angle, 
and not its amplitude. A knowledge of the singular stress fields near a three- 
dimensional re-entrant corner is, of course, required. 
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2.2 Fractal Dimension Evolution of Microcrack Net in Disordered 
Materials 


If the cumulative probability of the existence of microcracks whose 
lengths are longer than a specific value a, is 


P(a1) =[.pm&x (23) 


where p(x) is the probability density, and if we assume the form 


I 
è sara è sa 
p(x) = cNx +! = € a (pa) (24) 
N+2 
where C is a constant, and D= Ne ir then we obtain 
P(a,)=f  CNx “ax= [-x"] = Car". (25) 


With the cumulative distribution given by the above power law, the crack 
size distribution appears as a kind of fractal distribution with exponent N: 


POR InP(a;) — INC 


‘a 


N is only a measure of the degree of order in crack size distribution [6, 7]. 


(26) 


Alternatively, we can define D= [25] as the fractal dimension of 


N+1 
the crack size distribution in the specimen, so that the variation range of this 
fractal dimension is between 1 and 2, when N assumes positive values as 
shown in Fig. 3. Larger values of N correspond to more ordered situations, 
so that we can call N the “order parameter”. 

If the total number of cracks in the specimen is n; , and a; is any specific 
crack length, then 
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Nba, + Na, n n È (27) 
and 


no 2 fn Ste, fel (28) 


aa) 


where 
Ina, 


+ 
Ina, +InC-— Inn 


aa, 


IDIRÌ (29) 


According to the above definition, when D = 1, N tends to infinity and 
we have the ideal ordered material, just as when D = 1.5, N = 1 corresponds 
to the self-similarity defined by Carpinteri [6, 7], and the material turns out 
to be very disordered. 

If the total number of cracks n;, the minimum and maximum lengths of 
cracks, Amin and Amax , and the initial fractal dimension Dv are given as initial 


conditions, then from 
LE 1 1 
n n N ù La ’ (30) 


the constant C, corresponding to the given fractal crack size distribution, can 
be evaluated. 

For a given initial crack configuration, the actual crack size a; is obtained 
from the following estimation: 


1 1 i 3 


When i ranges from 1 to n, we get all crack lengths a; , 42, ..., dn. 


I 

For the sake of simplicity, uniform distributions are adopted for crack 
position and orientation. 

It is possible to study the evolution of the fractal  parameter 
corresponding to different degrees of disorder in crack size distribution 
during microcrack propagation. In the computer simulation, the total load is 
provided by different step increments. At each step, the number of cracks 
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with lengths greater than a given characteristic value 4; is computed, and the 
value of the fractal dimension is determined from equation (29). 

The evaluation of the total crack length during propagation is based on 
the assumption that the length of an initial crack with branching is 
effectively equal to the summation of the lengths of the branching cracks 
and that of the initial crack. 

From a theoretical point of view, the fractal dimension D should be the 
same for any reference crack length a; , if the distribution is exactly a fractal 
one. However, since we could only simulate a finite number of cracks and 
finite crack lengths, we take the average value over three evaluations for 
different crack lengths a, : 


PIE mn Amx + Inn cu +Anmn) 7 (32) 
i di ve 4 

In this way, we get three different fractal dimensions, of which we consider 

the average value for a specific loading level. 

In what follows, some examples are given for the purpose of studying the 
evolution of the fractal dimension of microcrack distribution during the 
propagation process, from the initial loading to the failure stage. These 
examples include different crack numbers and different initial fractal 
dimensions. In all cases the results show that the fractal dimension increases 
with the development of the microcrack net; i.e., the materials become more 
disordered with damage evolution. Although this analysis is only an initial 
attempt to study the evolution of the microcrack fractal dimension during 
the material damage process, the results are reasonable from the physical 
point of view since the longer cracks propagate more rapidly than the 
shorter, so that the microcrack size distribution becomes more dispersed. 

During crack propagation, some branching cracks may intersect with 
other primary or branching cracks. The assumption that the branching crack 
will stop propagating is adopted in the case of intersection. From the 
comparisons between different crack patterns obtained on the basis of arrest 
and non-arrest intersection crack assumptions, it is found that the results 
with arrest are more realistic. 

The first example deals with a square specimen of side 100 mm 
containing 200 microcracks with initial fractal dimension 1.40. The 
minimum initial crack length is 1 mm and the maximum is 10 mm. In each 
loading step, a stress increment of —4 N/mm? is applied. The original crack 
distribution and the subsequent crack patterns up to failure are shown in 
Figures 4.a to 4.f. The variation of the fractal dimension in terms of crack 
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size distribution is shown in Fig. 5, where the fractal dimension increases 
with the development of the microcrack net during the material damage 
process. This means that the material becomes more disordered when 
damaged. 

With the same number of cracks and the same initial minimum and 
maximum crack lengths, different initial fractal dimensions yield different 
crack size distributions in the specimen. The second example regards the 
evolution of the different initial fractal dimensions during the material 
damage process. The number of cracks is 100, and the initial minimum and 
maximum crack lengths are again 1 mm and 10 mm, respectively. Five 
initial fractal dimensions are considered, ranging from 1.263 to 1.423 (a 
stress increment of -4 N/mm? in each loading step is applied). Their 
evolutions are shown in Fig. 6. 


3. Incomplete Similarity and Application of Renormalization Group 
Theory 


Let us assume that the reacting section or ligament of a disordered 
material at peak stress could be represented as a fractal set of dimension 
2-d,, with O<d, < 1. The dimensional decrement d, may be due to 


the presence of voids and cracks and hence, generally, to a cross-sectional 
weakening. Let us consider two bodies, geometrically similar and made up 
of the same disordered material. If the ratio of geometrical similitude is 


equal to b and the renormalized tensile strength G,, is assumed to be a 


material constant and to have the physical dimensions [F] [L] (7%) , We 
have [8] 


(33) 


F, and F being the ultimate tensile forces acting on the two bodies, 


respectively. 
On the other hand, the apparent nominal tensile strengths are, 
respectively 


u tà ts Ou dr: (34) 
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where the latter, according to Eg. (33), becomes 
c,=0"b®*. (35) 


We can write the relationship between nominal strengths related to different 
sizes in logarithmic form 


Ino, = Ino! — d, Inb. (36) 


An alternative way to explain the decrease of the nominal tensile strength 
with specimen size is that of considering a sequence of scales of observation 
[9]. If the total force transmitted to the specimen is invariant with respect to 
the scale of observation, we have 


F=0;A;j=0,A,=..=G,An 


=0,An=0 STA 


A; == 0A, 


n+l 


where the first scale of observation could be the macroscopic one, with 
c,A,=6,A,A being the cross-sectional area, and the asymptotic scale of 
observation could be the microscopic one, with GA, = 0A , A' being 


the measure of the lacunar fractal set representing the damaged ligament. 
From the equality between the extreme members of Eq. (37), we get 


[A 
Ou “i u A .’ (38) 
and therefore 
bi È 
(07 i oi( b? ) (39) 


with b equal to the characteristic dimension of the cross section. From Eq. 
(39) we can get a generalization of Eq. (36): 


Ino, =lno,-d,Inb. (40) 
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Whereas Eq. (36) derives from the comparison between two specimens of 
different size-scale (simplified procedure), Eq. (40) derives from the 
comparison between two extreme scales of observation (cascade procedure 


or renormalization group). The equality of the intercepts provides 


* Bio è ; 
0,=0 @ , 1.e., the renormalized strength presents the same value as the 


nominal strength for the specimen of unit size. On the other hand, these two 
quantities present different physical dimensions and thus different units of 
measure. It is interesting to consider the graphical meaning of the 
Renormalization Group Procedure  (RGP), which is that of a 
counterclockwise rotation about the intercept (Fig. 7). 

Let us consider two bodies, geometrically similar and made of the same 
disordered material. If the ratio of geometrical similitude is equal to b and 


the renormalized fracture energy CI € is assumed to be a material constant 


and to have the physical dimensions [F] [L] (Hd) , we obtain [8] 








2xdg > (41) 


W, and W being the energies dissipated in the two bodies, respectively. 
On the other hand, the apparent nominal fracture energies are respectively 





E 1? ’ I) b? È 
where the latter, according to Eg. (41), becomes 
el AI dg 
pid. (43) 


We can write the relationship between fracture energies related to different 
sizes in logarithmic form: 


Inc@p = Inc@;! +d gInb . (44) 
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An alternative way of explaining the increase of fracture energy with 
specimen size is that of considering a sequence of scales of observation [9], 
as already shown for tensile strength. In an altogether analogous way, it is 
possible to find the following alternative expression: 


Incg = Inc95 +dgInb. (45) 


The comparison between Egs. (44) and (45) provides the equality 
Gp = GL between quantities with different physical dimensions. The 


graphical meaning of RGP is now that of a clockwise rotation about the 
intercept (Fig. 8). 

On the other hand, if the pioneering work of Griffith is revisited 
considering a fractal crack of projected length 2a, the fundamental relation 
of energy balance becomes [26]: 


2 

(O) 20 A +45 

Tama) = 2565 (0) DI 
and therefore 

o” UE des 

e da=95(1+dg)d “da , (47) 

From Eq. (47) we get 

ona” =(1+dg) GiE, (48) 


which represents the renormalized critical condition 
+\2 RE | 
(Ki ) vi (Kic) (49) 


The generalized stress-intensity factor K; presents the following physical 


dimensions: [F] [Lj (490)? . When dg = 0, we find again the classical 


relations of Griffith and Irwin. When dg =l,asa limit case, we find that 
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the stress singularity vanishes and Kj assumes the physical dimensions of 
stress. 


4. Scaling Laws for Tensile Strength and Fracture Energy of 
Heterogeneous Materials 


Both the dimensional decrement d, (tensile strength) and the 


dimensional increment d 4g (fracture energy) are experimentally always 


found to lie within the interval [0, 1/2 ]. At the macroscale, the size effects 


tend to vanish: d5 = dg =0, whereas, at the microscale, with the highest 


disorder, the size effects are the highest: d, =d g E 1/ 2. A material 


ligament of dimension 1.5 is equivalent to a Griffith crack distribution of 
self-similarity. On the other hand, a fracture surface of dimension 2.5 
(Brownian surface) is the least self-correlated locus (i.e., the most 
disordered), if we exclude overhangs. 

A generalization of the brittleness number may be provided by 
Dimensional Analysis [8-10] 


(50) 


If we take the reversal of the physical roles of strength and toughness to 
be absurd, the exponent of the characteristic size b must be positive: 


dara (51) 


This is a mesoscale inequality which must always be obeyed. We cannot 
exclude dg > 1/2 in the case of very intricate media or d, > 12 for very 
porous materials, provided that the inequality (51), however, holds. 

In physical reality, both material ligaments at peak stress and fracture 


surfaces after rupture can be considered as multifractals', of dimension 1.5 
and 2.5, respectively, at small scales, and dimension 2.0, at large scales. 


Some authors prefer the terms piecewise fractals or multiscale fractals. In any case, they should be 
defined in the framework of self-affinity instead of in that of self-similarity. 
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This means that, at large scales, the disorder is not visible, the size of the 
heterogeneities being limited. A transition from extreme disorder (slope 
12) to extreme order (zero slope) may therefore be evidenced in both 


diagrams, namely, Ino,, versus Inb, and Inc9€£ versus Inb (Fig. 9). 


On the basis of the foregoing arguments, two Multifractal Scaling Laws 
have been proposed for tensile strength and fracture energy, respectively 
[18-20, 27-29]: 


( 1/2 
O, = ri + a : (52.a) 





x l. =1/2 
GG (14 DÌ i (52.b) 


The asymptotic value of the nominal quantity (f, orC9£ ) is achieved 


only in the limit case of infinite sizes. 

A very convincing confirmation of Eq. (52.a) is given by the Brazilian 
tests carried out by Hasegawa e? al. [30]. Concrete cylinders, geometrically 
similar in two dimensions (with constant cylinder height of 500 mm), were 


tested in the scaling range 1:30 (Day =100mm, b,x = 3000 mm) . The 
Multifractal Scaling Law (MFSL) diagram is plotted directly in the 
bilogarithmic plane (Fig. 10), where the best fit from Bazant’s Size Effect 
Law (SEL) [31] is also shown for comparison. The computed best-fit values 
are: f,= 1.45 MPa and €,,=1992mm. Note that the asymptotic 


strength £, is 80% of the average ultimate tensile strength and only 56% of 


the strength of the smallest specimens. The correlation coefficient yielded 
by MFSL is R= 0.966, whereas SEL yields R= 0.663. In the case of 
tests characterized by wide scaling ranges, MFSL provides a better model of 
scaling behaviour than does SEL, since the concavity of the data distribution 
is clearly upwards. 

A confirmation of Eq. (52.b) has been provided by the three-point 
bending test results obtained by Elices et a/. [32]. 

The beam depth ranged from 50 to 300 mm (1:6). The nominal fracture 
energy was obtained from the total work of fracture, divided by the initial 
ligament. The application of MFSL is shown in Fig. 11. The best-fit values 


are respectively CH =1106 N/m and €, = 133 mm, The asymptotic 
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fracture energy is therefore almost 90% larger than the value for the smallest 
specimen (57 N/m). The correlation coefficient turns out to be R= 0.982. 
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Figure 1- Geometrically similar bodies with many random defects: constant defect 
size (a); defect size constant and proportional to the body size (b); defect size 
distribution of self-similarity. 
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Figure 2- Defect size distribution of self-similarity:  probability density (a); 
cumulative distribution (b). 
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Figure 3- Factal dimension D versus order parameter N. 
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Figure 4- Evolution of damage and increment of disorder in a compression test 
simulation. 
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Figure 5- Evolution of the fractal dimension during microcrack propagation. 
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Figure 6- Fractal dimension evolutions from different initial values. 
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Figure 7-Graphical interpretation of Renormalization Group Procedure for tensile 
strength. 
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Figure 8- Graphical interpretation of Renormalization Group Procedure for fracture 
energy. 
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Figure 9- Multifractal Scaling Laws (MFSL) for tensile strength (a) and fracture 
energy (b). 


196 ALBERTO CARPINTERI 


1.2 
1.0 
0.8 
0.6 
0.4 
0.2 
0.0 
-0.2 





Figure 10- MFSL versus SEL [30]. 





0 
O 50 100 150 200 250 300 350 400 


Figure 11- MFSL for best fit of fracture energy results [32]. 
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Aeroelastic Tailoring of Advanced 
Aircraft Wings Carrying External Stores 


Frank H. GERN', Liviu LIBRESCU” 


Abstract. /n this paper, structural and aeroelastic tailoring applied 
to advanced straight/swept aircraft wings carrying external stores is 
addressed. The wing structure is modeled as a laminated composite 
plate exhibiting flexibility in transverse shear and warping restraint 
effects. The equations of motion and boundary conditions are obtained 
via Hamilton's Variational Principle and application of generalized 
function theory. For a comprehensive representation of the stores, their 
static weights and inertia terms are considered. Modified strip theory 
aerodynamics is employed and the obtained eigenvalue/boundary value 
problems are being solved using the Extended Galerkin Method. The 
model is used to investigate the implications of external stores and ply 
angle variation upon divergence, free vibration, and flutter. Within the 
context of aeroelastic tailoring, it will be pointed out that already 
during the preliminary aircraft design phases, the influence of external 
stores attached to the wing structure has to be considered. 


Introduction 


Since modern civil and military aircraft wings are very often designed to 
carry external stores, this item has gained special importance in the 
consideration of wing aeroelasticity. In the case of transport aircraft these are 
mostly underwing carried stores such as large and heavy engines, fuel tanks, 
or tip stores like winglets. Experimental investigation of the winglet 
influence upon aeroelastic behavior revealed that most of the detrimental 
effects are due to the pure winglet mass and not due to aerodynamic 
interactions [1]. Store attachments to fighter aircraft mostly range within a 
wide variety of tank and missile configurations thus leading to literally 
hundreds of different store combinations. Very often, the stores exhibit 
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dramatic changes in mass and inertia properties during one single mission 
due to the launching of missiles or jettisoning of tanks [2]. 


Having in view that the tremendous variety of possible store 
configurations dramatically influences static and dynamic aeroelastic 
behavior of aircraft wings, the study of this problem has received 
prominence with the upcoming of variable sweep fighter aircraft only two 
decades ago [3], [4]. Since it is known that pylon mounted stores strongly 
influence dynamic wing characteristics, the store pitching modes are of 
special importance from an aeroelastic point of view [5]. With store pitching 
modes being in the proximity of the wings fundamental bending frequency, 
a critical aeroelastic coupling of the modes may occur, commonly referred to 
as wing-with-stores flutter. 


The possibilities of beneficially tailoring the aeroelastic properties of 
wing structures by the employment of composite materials have been 
pointed out by recognized aeroelasticians, (see e.g. [6], [7]), and were 
demonstrated by the development of the Grumman X-29 forward swept 
wing aircraft. But, as shown in Refs. [8] and [9], non-classical and often 
detrimental effects like warping inhibition and transverse shear flexibility of 
constituent materials are important parameters which have to be included in 
the aeroelastic model. 


As a result, for a successful design of the next generations of aerospace 
vehicles, and in order to eliminate 
the danger of the occurrence of any 
aeroelastic instability jeopardizing 
their imposed missions, it is 
imperative to have a perfect 
knowledge of the circumstances 
yielding the most 0 critical 
aeroelastic instability. Towards this 
end, as a necessary requirement, 
modeling of the flight vehicle 
should include the effect of 
distributed wing and tip stores, as 
well as consider the wing manufac- 
tured of composite materials, and 
investigate the flight in various 
speed regimes. 





line of aerodynamie center 





Figure 1- Geometry of the swept 
composite wing carrying external stores 
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Structural Modeling 


In order to investigate the effects of wing mounted stores upon 
aeroelasticity of advanced aircraft wings, a comprehensive structural model 
has been developed. Based upon the concept of a shear deformable plate- 
beam model, the wing structure is idealized as a laminated composite plate, 
each constituent layer featuring different ply angles, material, and thickness 
properties. The total number of the constituent layers is denoted by N. The 
reference plane of the composite structure is selected to coincide with the 


plane interface between the two contiguous layers r and r + l (1<r<N). 
Its points are referred to a Cartesian system of in-plane coordinates (xj, x,). 
The coordinate x} is normal to the plane (x, x,), its positive direction being 
the upward one (Figure 1). The xj and x, coordinates are referred to as 


chordwise and spanwise coordinates, respectively, whereas the reference 
plane is defined by x3 = 0. 


The structural model is developed postulating chordwise non- 
deformability as well as the following representation of the 3-D 
displacement components: 


U; (x, it} pensa) =; (x, ska 1) +xX34 (x, ME t), 
U, (x, saga t) =U (x, 308; t) +XY (x, Mai t), (la-c) 
erita aisi ha u3(x1,x3:1) 

Here (u,, u,, uz) denote the displacement components along the (x;, x), 
X3) coordinates and are associated with the points of the reference plane of 
the composite panel (x = 0); w, and y, denote the angles of rotation of a 
line element originally normal to the reference plane about the axes x, and 
x, respectively, while f denotes the time variable. Equations (1), well-known 
in the modeling of shear deformable plate and shell theories correspond to 
the first order transverse shear deformation theory (FSDT) (see e.g. Ref. [9]). 


In order to reduce the 3-D problem of the wing structure to an equivalent 
1-D one it is further postulated that 


Y, (x, 13055 1) = 9(x, ; t), 
Y, (x, FR 1) = f, (x, ;1) +x,8) (x, st), (2a-c) 
us(x, LOSE 1) = h(x, ; 1) _ (x, —%0 )e(x, 5 1) 
Equation (2c) describes the vertical displacement of a wing of rigid cross 
sections where h(x,:1) denotes the plunging displacement (positive upward) 
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of the wing cross-sections measured at the elastic axis, located at 
Xx = Xp (& xy (X,)). As a result, Egs. (1) become: 
U, = x30(x,;1), 
U,=u>(x,:1)+x;[f(03:0)+ argo (0451 JI (3a-c) 
U,= h(x,;1)-(x, 0 )e(x, ; 1) 


Consistent with Egs. (3), the strain components e;; are [10]: 
e,=U,,=0, 


Sergi 


=x,(0.+8g.) (4a-f) 


= 
(N°) 
i) 
3) 
Sn 
Il 
iS 
IS 
+ 
Sa eq 


Via E 2e3)=U,3 +U,,=0, 
Yx(& 2e,,)= SIA At a+ 40 +(x,0), 


In Egs. (4) and in the forthcoming developments, () j = 9()/0x;. 


For f, = -h3-(X9); and g, = 0), it results that y;3 = %3 = 0 which is 
consistent with the traditional assumption of an infinite stiffness of the wing 
structure in transverse shear (Kirchhoff°s Theory). 


The relevant equations of motion as well as the appropriate boundary 
conditions are obtained via Hamilton’s Variational Principle and application 
of generalized function theory in order to exactly consider the spanwise 
location and properties of the attached stores. To achieve a realistic 
representation of the store influence upon static and dynamic behavior of the 
entire system, static weight and inertia of the attached stores are modeled, 
yielding the energy functional of a composite wing carrying external stores: 


HI 
J=|(Uy —Ky+Ay-Ks+A5)dt (5) 

lo 
In equation (5) U denotes the strain energy, K represents the kinetic 
energy, and A the potential energy of body and surface forces. The indices W 
and S affecting the various quantities identify their affiliation to the wing and 
externally mounted stores, respectively. From the stationary condition à/ = 
0, consistent with equation (5) and adopting the Einstein summation 

convention, one obtains 
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dl =0= ja - [0,6U,,dt+{p(a,-U)6U dt+[o,sU, dal ‘9 
T t 5 


lo 


wing 


] 
+ fa. - x) [PATÙ 10 DO - MOD? + SUA dx, 


s_0 T, 





external stores 
In Eq. (6), dè denotes the variational operator, é,, represents Dirac’s 
distribution, whereas the superposed dots denote time derivatives. The terms 
underscored by a tilde represent prescribed quantities and the index s 
identifies the affiliation of the respective quantity to the external stores. # 
represents the component of the body force vector #/(per unit mass), whereas 
0;; denotes the respective component of the stress tensor. 


In the case of a straight wing, the pitching angle @ of the stores 
coincides with the wing pitching angle 0. By consideration of wings 
featuring the sweep angle A, the 3-D displacement quantities of the stores 
have to be referred to the wing coordinate system in order to express the 
stores’ main inertia properties which are in chordwise direction. This 
transformation yields 


Li 
Ue: 
Un (x - xo {@cos A — f, sinA — x,g, sin A] 


X-0cosfh—u, sinA—x, (f. +83 )sin A 


Il 


x30sinA+u,cosA+ x; (f. kan )cos A (7a-c) 


Employing the displacement and strain components from Egs. (3) and (4) 
and performing the indicated mathematical operations, the explicit form of 
equation (6) for the case of a composite aircraft wing carrying external stores is: 


0 


t 1 
d/=0= Jal {[Adw, + A,df, + A:dg, +(A- M ) 9 +(A+ L) dh]dx, + A; 


energy functional of the clean composite wing 
= Ym®8,(x, = x [(i. + NI, +x,t3 8, )òu, (8) 
= 2 (5) . 2 =; 2 29). “ 25)» fe 
+xi + +K,° sin NÈ, +[xx; +5 sin AK OsinAcosA-E" hsinA]d} 


PS 2 25) +. 2 pa O 25) + S CEI Si 
+x; [usi Ha? +K5° sin? A if Axe? +xK, sin A &-K 'OsinAcosA-E ÎrsinA)è, 
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- (- K°° f,sinAcosA-x,K}" è, sinAcosA+ (n? + K?" cos? Ab 

+ E hcos? A)60 È (- E° fi, sinA-x,E®&,sinA+E6cosA+ i)an] 

rr _=—=ell/P"®Òe:l!! W  NWMIPÒU”RU02ÒPE.,î‘” --—-}|;|)|.|e-.) ‘è ijiP jijijji iijfjfiblicE«cvcHYÒ:KKIODOIO:EDO]GO !‘ROROe 
kinetic energy of external stores 

ti Ys(x - x Mm ® sin Agf, +.x,M © sin Agg, — M ‘ cosA60 + 9h] 


aerodynamic forces acting on external stores 


+ im sò, \E sin Adf, +x,E}” sinAdg, — ES” cosAd0 + a] 
potential energy of external stores 


Herein, £,(9 denotes the distance between the center of gravity of the store 
and the wing elastic axis and Kl®) denotes the pitching radius of gyration of 


the store about its center of gravity, while g denotes the gravitational 
acceleration. The coefficients A, to Ag are displayed in Ref. [10] and, for the 


sake of completeness, are recorded here: 
A, =[r09 +F00 - 10%, - PASSA a ec 
A, = ia u Teo + FOO 1%, -1°9}, 212}, I 
PIC (74) - TO0 IO + FOO _]ODy, -1I9f, =19Ag;|, 
Ai = ba Lo si x Lp. Fa at e E (9a-f) 
(iS +199 -2x109 +19) +(109 400], 
A; = °° ta Foo - 109] Da FAI — x 799 6]. 
Age [7 8u, I, i (Tex. t a [7989], na [70°] 


l 
0,1 1,0) 0,0 
= [(76 Ta +xT5 ')so], 
In Equations (9), the one-dimensional generalized stress couples and the 
generalized body forces and mass terms, measured per spanwise unit length, 


are defined as: 


i Proggna E SV È O;;xX,"x3"dA (10) 


FE (o) H; MIE: 
N: l È X3 dA (11) 
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Collecting the terms associated with the respective variations du), dh, 
dg, 60, and dh, and having in view that these variations have to be arbitrary 


and independent, equation (8) yields the equations of motion as well as the 
boundary conditions. From the stationary condition d/ = 0 which concerns 
each instant belonging to the interval [#, ij], the equations of motion result 


as: 


du): A, +28(x 2-0 (+xÎ +10,8,)]=0 (12a-e) 
Halk Dil x in 1° (si 4(°+R2 sin? PA Ax +65 ind), 

— K° 6 sinAcosA — EVhsinA)+ M © sinA +mgE 0 sina]=0 
dg, A+Z5(% -#)x|- 


_k% sin A cos A — ESÎsinA)+ M  sinA+m® gE0® sina|=0 


a (x +x7 + K ® Sin? A) D- +? + 0° sin A) è 





60 : A+ Loox x" m°(- Ki "fi sinAcosA-x,K}" 8, sinAcosA 


+(47 x +K?"° cos? Alb +E® cos? A)- M ‘° cosA-m°° gE®® cosA| =0 
dh: A+g{x ate f.sinA-xE8, sinA+EdcosA+h) +1 +09 g| =0 


In addition, the geometrical and statical boundary conditions at the wing 
root and tip are obtained. For the case of a cantilevered wing, the boundary 
conditions at the root (x, = 0) are purely geometrical, expressed as: 


U,=U,, f1=f,: 8:=8,: 0=0, h=h (13a-e) 


In the case of a clean wing tip (no tip store), the boundary conditions at the 
wing tip (x, = /) are purely statical: 


. 0,0 0,0 
du: T,° = TO, 
IT (0,1) (0,1) 
d, Pi Ty = had 
dg,: ) den È: o sg (14a-e) 


+ mp(0,1) (1,0) _ (0,1) (1,0) 
60: To = {a =ITL sile 

. 7p(0,0) _ (0,0) 
oh: TI =TO 
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When considering a store located at the wing tip, the boundary conditions 
at x, = / change to kinetic ones thus taking into account the mass, inertia 


properties, and aerodynamics of the tip store: 

du, : TOO = -mP(î,+x, vo PUCE ) (15a-e) 
se E 

Tm xii +(x) +K° 


AT)A 


sin A)f, +(x,0;? + x KI0 sin Ali, 

6 sinAcosA-Elhsin Alt MsinA+m gEl” sin A 

dg,: di se - m°( (xs +(x) + KIT sin? A)f, + (110; +x 30 sinA)é, 
l'é sinAcosA— EDlisinA)+ Ml sinA+m gEl!® sin A| 

Spia Dr ele E cal =m"|-K FAL sinAcosA-x, Kik '&, sinAcosA 

n (x + K°"°° cos A) +E! hr cos a} M cosA-m°gEl" cosA 
h:19° = mP(- E® f,sinA-xE 8, sinA+E6cosA+h)+ LP +mP g 

Equations (12) through (15), considered in conjunction with Egs. (9) 
through (11), represent the equations governing the aeroelastic equilibrium 
of advanced composite aircraft wings laminated of anisotropic layers and 
exhibiting transverse shear flexibility and warping inhibition. In addition, 
they include arbitrarily distributed external stores in the wings spanwise and 
chordwise directions as well as the aerodynamics of the stores. These 
equations, expressed in terms of the unknown displacement quantities 
u(xg;t), fx; 1), QX; 1), 0 (x; 1) and h(x,; 1), yield a tenth order 
governing system of ordinary differential equations. 

Upon discarding the influence of the in-plane components of body forces 
F:°%, as well as in-plane rotatory inertia terms (i.e. terms /‘%%,, 10% fi, 
and /°" &,), the displacement quantity w(x,; 1) can be expressed in terms of 
fxa; 1), go(X5; 1), 0 (x; 1) and A(x,; 1) and thus be eliminated from the 


system. By this way, the system can be equivalently reduced to an eighth 
order differential equation system in terms of the unknowns f(x; 1), g2(X3; 


1), 8.(x,; 1) and h(x,; 1). 


Constitutive Equations 
The constitutive equations relating the generalized stress couples 7)" 


with the strain measures have been obtained in Ref. [10]. In order to 
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establish to study the effects of external stores within the context of 
aeroelastic tailoring, the composite wing is chosen to be fabricated of a finite 
number N of homogeneous layers. It is further assumed that the material of 
each layer is linearly elastic and that the bonding between the layers is 
perfect. The three-dimensional constitutive equations for a generally 
orthotropic elastic material can be expressed as 


Oi Q,, Q, di, 0 0 Qi € 
On, do O; d, OTTO, e E, 
Salo Gila 050 0A (16) 
03 | | 0 0 0 Q,, Og 0 Ya 
013 GA, 0, Oka Wa 
sal dx Dx DE 0 0° Via 


where Q, denote the transformed elastic coefficients associated with the 


k-th layer in the global coordinate system of the wing structure; Y;; = 2&;;, i#j 


and €;; denote the components of the strain tensor. 


Numerical Solution of the Static Aeroelastic System 
As previously stated, the Extended Galerkin Method has been applied for 


a numerical solution of the problem. The displacement field is represented as 
the sum of a finite number of mode shape functions as 


[1,8 (mM.0m.1M] =Y[F,,6,,7,,4,] n° (17) 
j=l 


where n denotes the non-dimensional spanwise coordinate (n= x,/l), 
while the constant factors /;, G;, T;, and 7; represent the modal amplitudes. 


Due to the structural complexity of the model and in order to establish 
trends, strip theory aerodynamics have been employed. However, it should 
be stressed that application of the Extended Galerkin Method allows a three- 
dimensional integration of the aerodynamic forces and moments along the 
wing span. Therefore, the sectional lift-curve slope and aerodynamic center 
remain arbitrary and thus may vary from section to section. 


For the static case, the aerodynamic terms / and / representing the 
sectional lift and aerodynamic torsional moment are expressed as [11]: 
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L(M) = In cA99,; 


18a,b 
M (N) FT q,ce4,8,; | 


where q,=pV, /2=qcos° A represents. the. dynamic pressure 
component normal to the leading edge, a, = 27AR/(AR+4cosA) denotes 
the lift curve slope coefficient corrected as to include the effects of a finite 
wing span and the wing sweep angle A. In equation (18), 0,y is the effective 
sectional angle of attack given by 


0, 


1 oh 
off Ah SR Ea (19) 


With @, denoting the angle of attack of the rigid wing structure, equation 
(19) including aerodynamic bending/twist coupling describes the well- 
known wash-in and wash-out effects of swept-forward and swept-back wing 
configurations, respectively. Employing Egs. (17), (18) and (19) into Eq. (8), 
an inhomogeneous matrix equation describing the static aeroelastic response 
of the wing is being obtained: 


[A]-[s]=[r] (20) 


The matrix A is composed of the (real) stiffness matrix of the wing but 
contains also the (real) aerodynamic influence quantities in terms of the 
dynamic pressure gy , whereas Ss represents the solution vector with 


s=[F---F,,G--G 
rates the inhnomogeneous part of the aerodynamic forces associated with the 
angle of attack of the rigid wing and the influence of the external stores. In 
order to obtain the divergence pressure, the determinant of matrix A has to 
be calculated, yielding a characteristic polynomial in g,. The smallest 
positive value of g, fulfilling the condition det A=0 represents the 


TT; HxH} The right-hand vector r_incorpo- 


n? 


divergence pressure (9, A i 

Since there is no contribution of the external stores to the matrix A, it 
becomes now evident that the influence of the stores on the divergence speed 
of the wing is immaterial. Experimental results obtained by Runyan and 
Watkins [12] reveal an identical behavior. 


Numerical Solution of the Dynamic Aeroelastic System 


For simulation of the dynamic aeroelastic system, the unknown functions 
are represented as 
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[10.8 0:0).0@)hm0T = [E] n e 1) 
j5l 


For the dynamic case, the aerodynamic terms / and /7 representing the 
sectional lift and aerodynamic torsional moment associated with an 
incompressible flow field are expressed as [11]: 


2,3) 1 oh 190 
Mp0 ha da ali far:Magr ae Luo} 


1 on ln 02) 
2,4)h 1 oh 190} 
M (N,1)= 7p@°b'4— My, +-=—Mgy;+bMy agi 
CREO I #*709n ron Mal 
In equation (22), Ly, Lypw» -.- » Meg represent the aerodynamic 


coefficients, whose expressions are displayed in Ref. [11]. 


For the Theodorsen function C(k), the approximation given in Ref. [13] 
was used, namely 


C(k) = F(K)+iG(k) 
_— 0.021573+0.210400% + 0.512607k? +0.500502%* (23) 
0.021508 + 0.251239% +1035378k? + k° 


-0.001995+ 0.327214% + 0.122397k? + 0.000146kÈ 
0.089318 + 0.934530% + 2.481481k? + £° 


with £ denoting the reduced frequency 


1 aa —] cereag vw 24) 
; V V cosA ( 


This representation finally leads to a complex eigenvalue problem 
expressed in matrix form as: 


[A]-@?[B]=0 (25) 


where A denotes the (real) stiffness matrix of the wing and B is the 
(complex) matrix representing the inertia terms of wing and external stores 
as well as the complex aerodynamic parameters. The real part of the 
complex valued quantity © represents the circular frequency of the 
oscillation, whereas its imaginary part constitutes the damping factor è. 


The implemented solution methodology is based upon the inversion of 
the complex matrix B and subsequent calculation of complex eigenvalues 
and eigenvectors of the obtained system matrix AB. The flutter speed is 
calculated in a fast converging iteration process rendering zero the imaginary 
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(damping) part of the complex eigenvalues. By this way, unlike in other 
models of similar structural complexity, the solution procedure neither 
involves a trial and error process [8], [14] nor requires to enforce 
convergence of the solution towards the flutter speed of the system [15]. 


Validation 


In order to verify the accuracy of the flutter analysis, a number of 
comparisons for several test cases were conducted. The first comparison was 
made with Goland’s cantilevered wing of AR = 6.67 as given in Ref. [16] 
and the subsequently appended correction of the flutter results in Ref. [17]. 
The results reveal that the predictions for flutter speed and flutter frequency 
provided by the present approach, namely 493.6km/h and 12.02Hz, are in 
excellent agreement with Goland’s exact results (494km/h and 11.25Hz). 


As another test case, the example of a straight aircraft wing of AR = 6.16 
with attached tip weights as given in Ref. [17] was investigated. In the 
calculations, two different chordwise positions of the center of gravity of the 
tip store were considered, namely £, = £,1%/c = 0.0 (case I) and £7 = 0.1 


(case II). It should be mentioned that the original flutter investigation by 
Goland and Luke also included the free body motion of the rigid fuselage, 
which is not subject of the present calculations. Nevertheless, the prediction 
of flutter speed and frequency of the first torsional mode for case I and case 
II is very accurate. 


In both cases, flutter occurs in the classical way as binary wing 
bending/torsion flutter due to frequency coalescence and vanishing 
aerodynamic damping of the first torsional mode at a flight speed of 1055 
km/h. In Refs. [15] and [18], it was observed that the first bending branch 
interacts with the rigid body mode and flutters at an even lower speed (994 
km/h). Since the free body motion of the fuselage was discarded in this 
investigation, there is no such interaction and, as a consequence, the first 
torsional mode yields the lowest flutter speed. 


In the present calculations, the warping inhibition of the cantilevered 
wing structure is taken into account. For this reason, flutter of the second 
bending mode which was observed by Housner and Stein [18] to be the most 
critical mode of instability (943 km/h) occurs only at very high air speeds 
[15]. The inclusion of the warping terms does not influence the stability 
behavior of the first bending and torsional modes (see, in this sense, also 
Ref. [19]). 


AEROELASTIC TAILORING OF ADVANCED AIRCRAFT WINGS... 213 


Numerical Illustrations 

In order to study the effects of external stores within the context of 
aeroelastic tailoring, the wing structure is considered to be manufactured of a 
graphite-epoxy composite material. In this context, the ply angle @ 
represents the counter-clockwise angle of rotation of the laminate with 
respect to the x,-axis of the wing. For a single layer composite wing, the 
material properties are given by: 


E;=30:10° psi; 


E,=0.75-10° psi ; 

G,,= 0.45-10° psi : 

G,3=0.37-10° psi ; 
v;,=0.25; 

p= 14.3-10° 1b sec?/in* 


For this type of 
material, the trans- 
verse shear flexibility 
parameter is 





0 20 40 60 80 100 120 140 160 180 


lulii de The effect of the 


ply angle @ upon the 
divergence speed is 
illustrated in Figure 2 
for a swept-back, 
swept-forward, and 
straight wing. The 
figure depicts the 
normalized diver- 
gence speed 


(4,),/ (4,), vs. the ply 


Figure 2- Normalized divergence speed vs. ply angle for 
swept and straight wings (R = 81, graphite epoxy; R + 0, 
Kirchhoff°s theory) 


1.05 


l 
(Gn)p.R=80 


(An)o.r-0 


0.95 
angle @, where 


(4.),=(4), se It can 
be seen that, with 


O 20 40 60 80 100 120 140 160 180 increasing forward 
Ply angle @ [deg] 





A=g= 





sweep, the range of 
ply angles for which 


Figure 3- Influence of the transverse shear flexibility vs. 
ply angle for different wing sweep angles 
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divergence is not 
critical, is decreased. 


Figure 2 and 
Figure 3 also high- 
light the effect of the 
transverse shear 
flexibility of the 
material. Despite the 
overall —detrimental 
effect of the 
transverse shear 
flexibility, Figure 3 
reveals that, for 
certain sweep angles 





0020 40 60 80 100 120 140 160 180 
Ply angle @ [deg] 


Figure 4- Normalized divergence speed vs. ply angle for 
swept and straight wings (FW = free warping, WR = 


A=0° | 
| A=-30° 


and 5 fiber warping restraint) 
orientations, consid- 
eration of the trans- 1.2 È er £ 









verse shear flexibility 
of the material may 
even result in an PE 
increase of the diver- rona 
gence speed. 0.8 | 





Figure 4 and 
Figure 5 display the 0.6 | 
influence of the 
warping restraint ef- 
fect, where 0.4 
1 (An)p.FW 020 40 60 80 100 120 140 160 180 
is the  divergence Ply angle @ [deg] 
speed for the free igure 5- Influence of the warping restraint effect vs. ply 
warping model and angle for different wing sweep angles 
(An)p.wr represents 
the divergence speed for the case in which the warping restraint effect is 
included in the simulation. 





The free vibration behavior of wings carrying external stores and 
exhibiting non-classical effects like transverse shear flexibility and warping 
inhibition can be obtained as a by-product of the dynamic analysis described 
in the above section. For this purpose, the airflow density is being set equal 
to zero, rendering the inertia matrix B real and therefore stating a self- 
adjoint eigensystem with purely real eigenvalues. 
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The influence of the warping restraint effect upon the eigenfrequencies of 
the clean wing is displayed in Figure 6. As it can be seen, neglecting the 
warping inhibition may lead to an underestimation of vibration frequencies, 
mainly within the vicinity of ply angles of around 90°. As expected, the free 
vibration behavior of the wing is symmetric with respect to a ply angle of 
90°of the single layer composite wing. 


Figure 7 illustrates the vibration behavior of a wing carrying different 
combinations of wing and tip stores exhibiting the following nondimensional 
properties: 





mi Tipstore) 
Hr Fo mi Vin® =0.1; 
mi Wingstore) 
Hy # mV sà 0.3 S 
(Tipstore) 
f e = =-04; 
[Hz] c 
(Wingstore) 
&= = 04; 
c 
x (Wingstore) 
n= = =05; 
0 20 40 60 80 100 120 140 160 180 l È 
Ply angle @ [deg] As depicted, vibra- 


Figure 6- Influence of warping restraint on the three tion frequencies are 
lowest eigenfrequencies for a straight wing (FW = free lowered due to an 
warping, WR = warping restraint) increase of the general- 
ized mass of the system 
passio by attaching external 

ni stores. 
TAR However, due to the 
chordwise offset of the 
center of gravity of the 
stores, the symmetry of 
the vibration behavior 
with respect to the ply 
angle is destroyed. As it 
can also be inferred 
O 20 40 60 80 100 120 140 160 180 fromFigure 7, the node 
i i lines of the higher order 
vibrational modes are 
shifted with respect to 
Ey=-0.4, my=0.5) changes in the  ply 


[Hz] 





Figure 7- Influence of external stores on the three lowest 
eigenfrequencies for a straight wing (u7=0.1, 4, = 0.3, er = 


216 


angle orientation as well 
by increasing the 
generalized mass [19]. 
The trend depicted in 
Figure 7 is consistent 
with the one obtained 
experimentally and 
theoretically by Lee [20]. 


For a non-zero airflow 
density, matrix B_ be- 
comes complex rendering 
the eigenvalues of the 
system to become com- 
plex quantities. Figure 8 
depicts the flutter speed 
of a clean wing vs. the 
ply angle @ for different 
sweep angles A. As it is 
shown, the maximum 
flutter speed of the wing 
is obtained for ply angles 
in the vicinity of 30° and 
is decreasing rapidly by 
approaching 0° as well as 
90°. This trend is in per- 
fect agreement with re- 


sults obtained by 
Housner and Stein using 
a finite-difference 


method [18]. 
To highlight the effect 
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Figure 8- Effect of ply angle @ and sweep angle A on 
flutter speed parameter Ag = Vp/b@y (influence of 


transverse shear flexibility) 
50 
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25 [4 H clean wing 
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Figure 9- Effect of ply angle @ on flutter speed parameter 
Àp = Vp/bay, (influence of external stores, 47 = 0.1, Ly = 


0.3, £e7= £y=-0.4, Ny=0.5) 


of the transverse shear flexibility on the flutter speed of the wing, Figure 8 
also depicts the flutter speed results obtained for a transverse rigid material 
(Kirchhoff theory, R+ 0). As it can be seen from Figure 8, neglecting the 
transverse shear flexibility leads to an overestimation of the flutter speeds for 


all ply angle orientations. 


The influence of different combinations of wing and tip stores is 
illustrated in Figure 9. Depending on the respective ply angle, the presence 
of external stores may lead to a decrease or increase of the flutter speed with 
respect to the clean wing. Moreover, as it can be seen from Figure 9, 
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changing the configuration of the external stores even leads to a change in 
the critical flutter mode for certain ply angles. For the structural 
configuration depicted in Figure 9, the attachment of a tip store renders 
mode 3 to become the more critical flutter mode for ply angles between 7° 
and 20°, whereas for all other configurations, mode 2 is the most critical 
flutter mode. 


Conclusions 


A well-encompassing structural model of aircraft wings composed of 
advanced composite anisotropic materials exhibiting transverse  shear 
flexibility and warping inhibition as well as incorporating arbitrarily 
distributed stores was developed. This model reveals its efficiency in 
approaching the static and dynamic aeroelasticity of complex wing/store 
configurations. Results obtained for Goland’s and for Goland and Luke”s 
wings via Extended Galerkin’s Method show very good agreement with 
solutions given by other authors and obtained via completely different 
approaches to the problem, such as analytical solution methodologies and 
finite difference. The capability of the model to accurately predict flutter and 
free vibration behavior of wings carrying external stores is underlined by 
comparison with theoretical and experimental results for standard unswept 
metallic wings and composite wings, respectively. 


The assumption of a non-shear deformable wing structure (R = E/G’ > 
0) overestimates divergence and flutter speeds as well as flutter frequencies. 
For this reason, in order to obtain reliable results for static and dynamic 
aeroelasticity of advanced composite aircraft wings, transverse shear 
flexibility has always to be taken into account. Since wing divergence is only 
influenced by wing stiffness parameters, the effect of external stores upon 
the divergence speed has been shown to become immaterial. 


Depending on the wing/store configuration, the implications of external 
stores on flutter speed of aeroelastically tailored composite wings can be 
very complex, in the sense that they can be of beneficial or detrimental 
nature. In addition, it has been shown that different store combinations may 
also lead to different critical flutter modes, depending on the ply angle 
orientation of the composite material as well. 


Ongoing activities are intended to evaluate the implications of the free 
body motion of the fuselage which in some instances may influence 
aeroelastic instability of forward-swept wing configurations by coupling of 
the first symmetric bending mode with the aircraft’s short period mode. This 
issue is envisaged to be approached in conjunction with the effects played by 
wing mounted stores and tip masses. Inclusion of the free body motion in the 
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aeroelastic study will require consideration of the mass and inertia terms of 
the fuselage resulting in a change of the geometrical boundary conditions at 
the wing root into kinetic ones. Consideration of the free body motion will 
enable one to approach also the asymmetrical flutter induced by non- 
symmetry of the store configurations on both sides of the wings. 


In the present analysis, the external stores are regarded to be rigidly 
attached to the wing structure. Therefore, the consideration of pylon 
flexibility will be an important step towards a more realistic representation 
of store influence upon wing aeroelasticity. Furthermore, in order to account 
to the structural complexity of the model, more refined aerodynamics 
(Doublet Lattice) are intended to be incorporated for upcoming analyses. 
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